How to find all neighbours of an element in N-dimensional matrix

I just can't wrap my head around how to find all neighbours of an element in a N-dimensional matrix. After reading through the documentation on spatial searching I think it could be done using the Delaunay triangulation. However, these topics are still a bit too advanced for me, so any help would be appreciated.

[dave]

Update: Thanks everybody for your suggestions. After all Shashank's hint with the knnsearch turned out to be the easiest solution in my case, since I have the statistics toolbox. If you don't have it, the solutions provided by Matt J and Image Analyst are also very useful.

Answer 1

[Matt J]

What I often like to do is make an offset mask, as follows

    s=size(yourmatrix);
    N=length(s);

    [c1{1:N}]=ndgrid(1:3);
    c2(1:N)={2};

    offsets=sub2ind(s,c1{:}) - sub2ind(s,c2{:})

Now, for any linear index in your matrix L, you can get all its neighbors by doing

    neighbors = yourmatrix(L+offsets)

It won't work at the edges of the matrix, but you can take care of that by padding the edges first.

[dave]

Thanks for the suggestion - however, I'm getting a "Subscript indices must either be real positive integers or logicals." error for the last line and I'm not sure why..

[Matt J]

It's because you are applying it to a point at the edge of the array, where you have fewer neighbors than in the array's interior. See my remarks above about padding as a possible solution for those points.

[dave]

So now I've padded the matrix with a padSize of 100 along all dimensions and the error's still there...

[Matt J]

What is the size of your array and what is the index value you are using for L?

[Matt J]

I think you need an example to see what this method is doing. Suppose yourmatrix is

yourmatrix =

       1     6    11    16
       2     7    12    17
       3     8    13    18
       4     9    14    19
       5    10    15    20

When you run my code on this small example, you should obtain

 offsets =

    -6    -1     4
    -5     0     5
    -4     1     6

So, now, if I want the neighborhood centered at the "14", I would do

    >> yourmatrix(14+offsets)

    ans =

         8    13    18
         9    14    19
        10    15    20

If, however, I do

>> yourmatrix(1+offsets)
Subscript indices must either be real positive integers or logicals.

it cannot work because the "1" is not in the interior and doesn't have a full 3x3 neighborhood. As a solution, I create a "picture frame" around the original image,

yourmatrix =

       0     0     0     0     0     0
       0     1     6    11    16     0
       0     2     7    12    17     0
       0     3     8    13    18     0
       0     4     9    14    19     0
       0     5    10    15    20     0
       0     0     0     0     0     0

Now, when I recompute "offsets", I obtain

offsets =

      -8    -1     6
      -7     0     7
      -6     1     8

And if I want the neighborhood surrounding the "1", I just do

 >> yourmatrix(9+offsets)

ans =

     0     0     0
     0     1     6
     0     2     7

Tada!!

[dave]

I see, thanks...

[Sotirios]

Matt J, I have a question if you can help me. I have done what you explained so far and I have created a 3D matrix containing in each index the neighbors of the initial matrix. Now I would like to manipulate this matrix (i.e. I would like to sum the values of each element that equals to each neighbor). However, since the matrix I have created also contains the picture frame with zeros it's hard for me to scan this matrix. What I mean its easy to set a start point in the loop (i.e N+2) but it's difficult to avoid the elements that belong to the last column of my matrix. Do you might have any suggestions? Thanks in advance.

[Matt J]

Hi Sotirios,

I don't completely understand the operation you are trying to do. If the question is how to avoid neighborhoods not completely contained in the image, one solution might be to not avoid them at all. Set the picture frame to Inf and perform the operation on all neighborhoods. Results that are Nan or Inf, you then throw away.

Another solution would be to only generate neighborhoods for interior pixels

 T=true(size(your3Darray));
 T([1,end],[1,end],[1,end])=0; %exclude edges
 L=find(T);

[Sotirios]

Yes that was my question. The operation has to do with the Ising model, where it is assumed that there is strong correlation between neighbor pixels in an image. But you answer to my question, I haven't thought of setting the picture frame to Nan.

Thanks again!

[Matt J]

I haven't thought of setting the picture frame to Nan.

Just bear in mind that NaN processing can be slower than Inf processing. I suggested setting the picture frame to Inf instead of NaN for that reason.

[Sotirios]

Alright, thank you for pointing it out because I didn't know it makes any difference.

[Roisin Loughnane]

This is a really helpful answer, and an informative use of offsets.

I am trying to use this in a loop over a large matrix containing NaN elements.

The problem is that the values of L are linear indices and when I try to loop over these I get the error:

Subscript indices must either be real positive integers or logicals.

How do I work around this?

[Matt J]

Does this happen even when you pad the array, as described above?

[Image Analyst]

Roisin, the FAQ explains the error well: http://matlab.wikia.com/wiki/FAQ#How_do_I_fix_the_error_.22Subscript_indices_must_either_be_real_positive_integers_or_logicals..22.3F

Answer 2

[Shashank Prasanna]

If you have the Statistics Toolbox installed, there are couple of nearest neighbor searching tools that you might find very useful. Although some of these concepts may be advanced for a casual user, the functionality itself is very quite easy to use, that's really the power of MATLAB:

One way to speed up higher dimensional neighbor searching is to use an optimized data structure such as k dimensional trees, and this is easy to do so as below:

kdtreeNS = KDTreeSearcher(x);
[n,d]=knnsearch(kdtreeNS,x);

Here is the page from the doc:

http://www.mathworks.com/help/stats/kdtreesearcher.knnsearch.html

[dave]

Thanks for the suggestion...Yes, I do have the Statistics Toolbox, but after reading through the doc for a while I'm still not sure how to apply KDTreeSearcher() to my n-dimensional matrix, since it obviously only accepts 2-dimensional input (X). What am I doing wrong?

[Shashank Prasanna]

Michael, I misunderstood your question. I made the assumption that you were referring to the columns as dimensions in |R^n space .

None of these function work on multidimensional matrices. If you have data in some n dimensional space each instance or observation can be represented as a vector with n entries. I am not sure how your data is represented but if you are able to represent it in a why described you can use all the functionality that I mentioned to do your neighbor search.

Answer 3

[Image Analyst]

You just take the index, the index plus one, and the index minus one, for every other dimension, but exclude the index of where you're at. For example if you have a 3D matrix, and you're at x=3, y=6, and z=9, you'd have all permutations of x in [2,3,4] with y in [5,6,7] with z in [8,9,10] but don't include the point itself with x=3,y=6, z=9. So that's 3*3*3-1 neighbors in 3D. For N dimensions you'd have 3 * 3 * 3 * ....(a total N times) * 3 -1 neighbors. So 2D gives 3*3-1 = 8 neighbors, 3D gives 26 neighbors, 4D gives 80 neighbors, etc.

[Image Analyst]

By the way, you mgiht like this video. http://www.youtube.com/watch?v=JkxieS-6WuA. It does a decent job of explaining all the dimensions up to the tenth dimension.

[dave]

I see.. I might be wrong, but doesn't this solution get quite exhaustive once the dimension gets larger? Can't it be done in a more efficient way that uses one of the built-in spatial searching functions?

[Image Analyst]

Yes it does. I believe though that you can carve out the box using indexing, for example:

subHyperVolume = hyperVolume(4:6, 18:20, 9:11, 7:9, 1:3);

though this sub volume would have the "center" element.

[Jan]

When I speak of a larger number of dimensions, I mean thousands to millions. Then interesting effects occur, e.g. that almost all points have almost the same distance. And the ratio between the volume of a hyper-cube and the embedded (or is it called "included") hyper-sphere has a surprising relation to the number of dimensions also.

But such large spaces are very hard to implement in Matlab, because the number of indexed elements explodes.

[dave]

@ImageAnalyst: I like your suggestion with the (sub-)hypervolume. Here's what I have so far:

indices = cell(1, length(size(hyperVolume)));      
for i = 1:length(hyperVolume(:)),                 
    [indices{:}] = ind2sub(size(hyperVolume), i);
     %...           
end

This way I iterate over every element of the hypervolume and store the indices of the current element. The only thing that's still unclear to me is how to get from the cell containing the indices to your line which extracts the subvolume:

endsubHyperVolume = hyperVolume(4:6, 18:20, 9:11, 7:9, 1:3);

Any idea?

[Image Analyst]

I'd use loops over each dimension, so for 5 dimensions like my example

for d1=1:d1max
for d2=1:d2max
for d3=1:d3max
for d4=1:d4max
for d5=1:d5max
subHyperVolume = hyperVolume(...
   d1-1:d1+1,...
   d2-1:d2+1,...
   d3-1:d3+1,...
   d4-1:d4+1,...
   d5-1:d5+1,...
);
end
end
end
end
end

Answer 4

[Sean de Wolski]

No need to reinvent the wheel:

http://www.mathworks.com/matlabcentral/fileexchange/29330-neighbour-points-in-a-matrix