Calculate shortest distance to selected points in 3D matrix in matlab

Question

I have a very large 3D image stored in a matrix (approx. 500x500x40 voxels) in Matlab. In this matrix about 30000 points are selected (suppose randomly). Suppose the selected voxels have a value of one and the non-selected points are zero. Now I need to calculate for every voxel in the entire 3D-image the Euclidian distance to the nearest selected point.

So for example in 2D, given a 4x4matrix:

    selection = 0 0 1 0
                1 0 0 0
                0 0 0 1
                0 0 0 0

The corresponding distance matrix would be:

    distance = 1  1  0  1
               0  1  1  1
               1  √2 1  0
               2  √5 √2 1

Is there an efficient way to do this, in terms of both time and memory?

How do you define your distance? Euclidian? Are the 3 dimensions of your matrix representing coordinates in 3 directions? — BillBokeey
Distance is indeed Euclidian. The 3D matrix is a medical scan, so every "point" in the matrix is in fact a voxel. — sdbonte
There are many possibilities, see i.e.: en.wikipedia.org/wiki/Nearest_neighbor_search . I was thinking about octree's, but I have no relevant experience with it. — Bernhard
Reshape your matrix into a 250000 x 40 matrix, then sample from this matrix to get a 30000 x 40 matrix of points, then use the above marked duplicate. — rayryeng
@sdbonte Right it didn't hit me about the memory storage until you mentioned it and I apologize for not seeing it properly before flagging. I'll remove the duplicate. I also thought that each 3D slice at one spatial coordinate was a data point, not the actual voxel itself. Because of the huge memory requirement, you may have to resort to storing on disk or splitting up the operation over multiple processors / machines. Either way, it's going to be rather slow. — rayryeng

Hans Hans · Accepted Answer · 2016-06-06T12:22:15

If you don't need to know which combinations had which distances, then you could calculate the 3D cross-correlation.

To illustrate this in 2D, take the following matrix and calculate the 2D correlation as described in the reference

M =
 0     0     0
 1     0     0
 1     0     0

convn(v,v(end:-1:1,end:-1:1)) =
 0     0     0     0     0
 0     0     1     0     0
 0     0     2     0     0
 0     0     1     0     0
 0     0     0     0     0

One can read off the distances between the ones because the correlation matrix here can be understood as the differences in indices. The central column in convn means the horizontal distance is zero. Likewise, the middle row means zero vertical distance. As a result, the central value gives you the amount of zero distance values, which is the sum of ones in your matrix. The two ones correspond to a vertical distance of 1 between the ones in M. One combination has a positive distance and the other combination has a negative distance.

As a result, you now have all the distances, but they contain the horizontal and vertical directions too. But you can still process that as you wish.

A way to calculate all possible squared distances without a for loop would be

n = size(M,1) = 
 3

tmp = repmat([-(n-1):(n-1)].^2,2*n-1,1)
d2 = tmp+tmp' =
 8     5     4     5     8
 5     2     1     2     5
 4     1     0     1     4
 5     2     1     2     5
 8     5     4     5     8

Both matrices together contain basically the histogram of distances.

Calculate shortest distance to selected points in 3D matrix in matlab

2 Answers