sdbonte
sdbonte

Reputation: 102

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

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?

Upvotes: 1

Views: 337

Answers (2)

Hans
Hans

Reputation: 46

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.

Upvotes: 1

Hans
Hans

Reputation: 46

If your points are given as coordinates X = n x 3, you can use

D = pdist(X,'euclidean')

to efficiently calculate all combinations of distances.

Upvotes: 0

Related Questions