Avoiding sub2ind and ind2sub

Question

I need to access several indices around a certain point in 3D. For example, for point (x1,y1,z1) I need to get all the indices of its 3x3x3 neighborhood such that (x1,y1,z1) is centered. For neighborhood of size 3, I do it with

 [x,y,z] = meshgrid(-1:1,-1:1,-1:1);
 x_neighbors = bsxfun(@plus,x,x1);
 y_neighbors = bsxfun(@plus,y,y1);
 z_neighbors = bsxfun(@plus,z,z1);

Here, I center x1,y1,z1 to (0,0,0) by adding the distances from (x1,y1,z1) to any point in the 3x3x3 box.

that gives me the coordinates of (x1,y1,z1) 3x3x3 neighborhood. I then need to turn them into linear indices so I can access them:

 lin_ind = sub2ind(size(volume),y_neighbors,x_neighbors,z_neighbors);

that is costly in what I do.

My question is, how to avoid sub2ind. If inx is the linear index of (x1,y1,z1),

inx = sub2ind(size(volume),y1,x1,z1);

how can I find the 3x3x3 neighborhood of the linear index by adding or subtracting or any other simple operation of inx?

Answer 1

As long as you know the dimensions of your 3D array, you can compute the linear offsets of all the elements of the 3x3x3 neighborhood. To illustrate this, consider a 2D example of a 4x5 matrix. The linear indices look like this:

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

The 3x3 neighborhood of 10 is [5 6 7 9 10 11 13 14 15]. The 3x3 neighborhood of 15 is [10 11 12 14 15 16 18 19 20]. If we subtract off the index of the central element, in both cases we get [-5 -4 -3 -1 0 1 3 4 5]. More generally, for MxN matrix we will have [-M-1 -M -M+1 -1 0 1 M-1 M M+1], or [(-M+[-1 0 1]) -1 0 1 (M+[-1 0 1])].

Generalizing to three dimensions, if the array is MxNxP, the linear index offsets from the central element will be [(-M*N+[-M-1 -M -M+1 -1 0 1 M-1 M M+1]) [-M-1 -M -M+1 -1 0 1 M-1 M M+1] (M*N+[-M-1 -M -M+1 -1 0 1 M-1 M M+1])]. You can reshape this to 3x3x3 if you wish.

Note that this sort of indexing doesn't deal well with edges; if you want to find the neighbors of an element on the edge of the array you should probably pad the array on all sides first (thereby changing M, N, and P).

Answer 2

Just adding the (generalized) code to @nhowe answer: This is an example for neighborhood of size 5X5X5, therefore r (the radius) is 2:

ns = 5;
r = 2;

[M,N,D] = size(vol);
rs = (1:ns)-(r+1);
% 2d generic coordinates:
neigh2d = bsxfun(@plus, M*rs,rs');
% 3d generic coordinates:
pages = (M*N)*rs;
pages = reshape(pages,1,1,length(pages));

neigh3d = bsxfun(@plus,neigh2d,pages);

to get any neighborhood of any linear index of vol, just add the linear index to neigh3d:

new_neigh = bxsfun(@plus,neigh3d, lin_index);