Voxel neighborhood indexing - Detecting "out of bounds" in 26 neighbor access with linear indexing

Question

Well, I don't know how to describe my problem with a title, I hope the one I got is correct.

I have a matrix (Min the example below) that is a 3D image, composed, in this case, by 11x11x11 voxels (I made it logical just for easiness, and size is just an example also).

In my code, I need to reach the 26 neighbors of some voxels, and for that I use some fancy linear indexing found in: http://www.mathworks.com/matlabcentral/answers/86900-how-to-find-all-neighbours-of-an-element-in-n-dimensional-matrix

The problem is that if the point is in the "boundary" of M some out of bounds values are tried to be accessed, and that will generate an error.

To solve this problem, a good approach would be to create a boundary around M making it +2 size in every dimension, and populate that with zeros, however I really would like to avoid changing M, as my code is quite more complex that the one in the example.

I cant find any way of doing it, i'm a bit stuck here. Any suggestion?

EDIT: @Dan answer works, however I would like to see if there is a possible solution using this linear indexing method.

% Example data
M=round(randn(11,11,11))~=0;

% Fancy way of storing 26 neigh indices for good accesing 
s=size(M);
N=length(s);
[c1{1:N}]=ndgrid(1:3);
c2(1:N)={2};
neigh26=sub2ind(s,c1{:}) - sub2ind(s,c2{:});

point=[5 1 6];

% This will work unless the point is in the boundary (like in this example)
neighbours=M(sub2ind(s,point(1),point(2),point(3))+neigh26) 

Answer 1

Is that linear indexing stuff essential? Because it's pretty easy to handle boundary conditions is you use subscript indexing and min and max like this:

p = [5, 1, 6];

neighbourhood = M(max(1,p(1)-1)):min(p(1)+1,end),
                  max(1,p(2)-1)):min(p(2)+1,end),
                  max(1,p(3)-1)):min(p(3)+1,end))

%// Get rid of the point it self (i.e. the center)
neighbours = neighbourhood([1:13, 15:end])

This way you can also easily generalize this if you want a broader neighbourhood:

p = [5, 1, 6];
n = 2;
neighbourhood = M(max(1,p(1)-n)):min(p(1)+n,end),
                  max(1,p(2)-n)):min(p(2)+n,end),
                  max(1,p(3)-n)):min(p(3)+n,end))

%// Get rid of the point it self (i.e. the center)
mid = ceil(numel(neigbourhood)/2);
neighbours = neighbourhood([1:mid-2, mid+1:end])

or if you liked to keep the cube shape then maybe:

neighbours = neighbourhood;
neighbours(mid) = NaN;

If you want to use this many times in your code it's probably best to refactor it as an m-file function that just returns the indices:

function ind = getNeighbours(M,p,n)
    M = zeros(size(M));
    M(max(1,p(1)-n)):min(p(1)+n,end), max(1,p(2)-n)):min(p(2)+n,end), max(1,p(3)-n)):min(p(3)+n,end)) = 1;
    M(p(1), p(2), p(3)) = 0;
    ind = find(M);
end