0

Ok so I've got two large 3D binary arrays. I want to calculate the minimum distance between the surface points of these two structures to get an output of the distances between the perimeter voxels of the objects displayed in the array data.

I know I'll have to use the bwperim() function; i.e. perim_a = bwperim(a,6); & perim_b = bwperim(b,6); and then feel like I should be able to use the norm() function to do something like norm(perim_a - perim_b); however I just keep getting the error 'Input must be 2d'. Is there a way to apply the norm() function to 3D data?

user3168953
  • 153
  • 2
  • 10

1 Answers1

2

If I understand correctly, you have two arrays of voxels in 3D that represent the perimeter - and you want the distance between voxels of one and the other. My approach would be something like this (assuming A and B are the initial volumes, and that they are the same size):

Aperi = bwperim(A,6);
Bperi = bwperim(B,6);

Now we need the coordinates of the non-zero voxels:

sz = size(A);
[xx yy zz] = meshgrid(1:sz(1), 1:sz(2), 1:sz(3));

indxA = find(Aperi==1);
posA = [xx(indxA(:)) yy(indxA(:)) zz(indxA(:))]; % a Nx3 matrix of x,y,z positions
indxB = find(Bperi==1);
posB = [xx(indxB(:)) yy(indxB(:)) zz(indxB(:))]; % a Mx3 matrix - note N ~= M

Now you have two 2D matrices and you can take the difference between points - but you don't know which pairs you want to compare so you need a bit more work:

delta = bsxfun(@minus, reshape(posA, 1, [], 3), reshape(posB, [], 1, 3) );

Now delta has dimensions [M N 3], with differences between each pair of points in posA and posB. The distance is

distance = sqrt(sum(delta.*delta, 3));

And for each point in A, the distance to the closest point in B is

closest = min(distance);

which will be a [1 N] row vector.

To put these values back into the original matrix, you need to go back to the indices you had in the beginning:

distMatrix = zeros(size(A));
distMatrix(indxA) = closest;

And I think you will now have the distance of the closest point on the perimeter of B to each point on the perimeter of A mapped to the surface of A.

Let me know if this gets you anywhere to the answer you were looking for...

Floris
  • 45,857
  • 6
  • 70
  • 122
  • I keep getting this error: `Error using reshape Product of known dimensions, 3, not divisible into total number of elements, 41418752. Error in D_surf_WP (line 47) delta = bsxfun(@minus, reshape(struc_a, 1, [], 3), reshape(struc_b, [], 1, 3) );` Any ideas how to solve? – user3168953 Feb 01 '14 at 14:11
  • also I'm nor sure I understand why you create the `'posA/b'` variables – user3168953 Feb 01 '14 at 15:49
  • I think the `reshape` should be on `posA` and `posB` - sorry on iPhone here can't check properly. Will take closer look later if that doesn't solve it. – Floris Feb 01 '14 at 17:35
  • @user3168953 - I see you accepted the answer. I assume that means my edit solved your two problems, and that it was my error (A instead of posA). As such I have now edited the answer. Please let me know if that was indeed correct. – Floris Feb 01 '14 at 22:00