Get the neighbors of a matrix element

Question

I have a matrix and for each element I want to get the index of its surrounding elements. All these results have to be stored into a matrix in the following way. Each row of the matrix corresponds to a matrix element and each of the columns of this matrix contain s the neighbor indexes. For example, for a 4x4 matrix we will get a 16x8 result array. Some of the matrix elements do not have 8 neighbors.

There is an example, I think it is working, I there any way to avoid for loop?:

ElementNeighbors = [];
for n = 1:numel(Matrix)
    NeighborsMask = [ n-1 n+1 n+size(matrix,1) n-size(Matrix,1) n-size(Matrix,1)-1 n-size(Matrix,1)+1 ...
        n+size(Matrix,1)-1 n+size(Matrix,1)+1 ];

    ElementNeighbors = [ElementNeighbors ; NeighborsMask ];
end
ElementNeighbors (ElementNeighbors ==0|ElementNeighbors <0) = NaN;

I don't understand the description of the desired output. Can you give an example line? — Daniel
There are solutions for this problem all over (e.g., here). What about the elements on the edge of the matrix that have fewer than eight neighbors? — horchler
That I was thinking about, a mask. For the border elements there will not be 8 neighbors. — Thoth

Floris Floris · Accepted Answer · 2014-01-30T18:21:05

Given the linear indices of a matrix M(n,m), you can convince yourself that the top left neighbor of element M(i,j) = M(i-1, j-1) = M(i-1 + n * (j-2))

In "linear index" space that means the offset of this element is

-n-1

Doing this for all other locations, we find

-n-1 | -1 | n-1
-n   |  x | n    => [-n-1, -n, -n+1, -1, +1, +n-1, +n, +n+1]
-n+1 | +1 | n+1

Thus you can create a vector offset with the above values (replacing n with the first dimension). For example, if M is (5x4), then

offset = [-6 -5 -4 -1 1 4 5 6];

You then create all the indices:

indices = bsxfun(@plus, (1:m*n), offset(:));

bsxfun is a cool shorthand for "do this function on these elements; where one element has a singleton dimension and the other doesn't, expand accordingly". You could do the same with repmat, but that creates unnecessary intermediate matrices (which can sometimes be very large).

That command will create a (8 x m*n) matrix of indices of all 8 neighbors, including ones that may not really be the neighbors... something you need to fix.

Several possible approaches:

pad the matrix before you start
don't care about wrapping, and just get rid of the elements that fall off the edge
create a mask for all the ones that are "off the edge".

I prefer the latter. "Off the edge" means:

going up in the top row
going left in the left column
going down in the bottom row
going right in the right column

In each of these four cases there are 3 indices that are 'invalid'. Their position in the above matrix can be determined as follows:

mask = zeros(size(M));
mask(:,1) = 1;
left = find(mask == 1);
mask(:,end) = 2;
right = find(mask == 2);
mask(1,:) = 3;
top = find(mask == 3);
mask(end,:) = 4;
bottom = find(mask == 4);

edgeMask = ones(8,m*n);
edgeMask(1:3, top) = 0;
edgeMask([1 4 6], left) = 0;
edgeMask([3 5 8], right) = 0;
edgeMask(6:8, bottom) = 0;

Now you have everything you need - all the indices, and the "invalid" ones. Without loops.

If you were feeling ambitious you could turn this into a cell array but it will be slower than using the full array + mask. For example if you want to find the average of all the neighbors of a value, you can do

meanNeighbor = reshape(sum(M(indices).*edgeMask, 1)./sum(edgeMask, 1), size(M));

EDIT re-reading your question I see you wanted a M*N, 8 dimension. My code is transposed. I'm sure you can figure out how to adapt it...

ATTRIBUTION @Tin helpfully suggested many great edits to the above post, but they were rejected in the review process. I cannot totally undo that injustice - but would like to record my thanks here.

EXTENDING TO DIFFERENT REGIONS AND MULTIPLE DIMENSIONS

If you have an N-dimensional image matrix M, you could find the neighbors as follows:

temp = zeros(size(M));
temp(1:3,1:3,1:3) = 1;
temp(2,2,2) = 2;
offsets = find(temp==1) - find(temp==2);

If you want a region that is a certain radius in size, you could do

sz = size(M);
[xx yy zz] = meshgrid(1:sz(1), 1:sz(2), 1:sz(3));
center = round(sz/2);
rr = sqrt((xx - center(1)).^2 + (yy - center(2)).^2 + (zz - center(3)).^2);
offsets = find(rr < radius) - find(rr < 0.001);

You can probably figure out how to deal with the problem of edges along the lines shown earlier for the 2D case.

Untested - please see if you notice any problems with the above.

Get the neighbors of a matrix element

1 Answers