Индекс индекса Voxel - обнаружение «вне границ» в 26 соседнем доступе с линейной индексацией

Question

Ну, я не знаю, как описать мою проблему с заголовком, я надеюсь, что тот, который я получил, верен.

У меня есть матрица (M в приведенном ниже примере), который представляет собой изображение 3D, состоит, в данном случае, по 11x11x11 вокселей (я сделал это логично только для простоты, и размер это всего лишь пример и).

В моем коде, я должен достичь 26 соседей некоторых вокселей, и для этого я использую некоторые фантазии линейной индексации найти в: http://www.mathworks.com/matlabcentral/answers/86900-how-to-find-all-neighbours-of-an-element-in-n-dimensional-matrix

Проблема в том, что если point находится в «границе» из M Пробуется доступ к некоторым значениям вне границ, и это приведет к возникновению ошибки.

Чтобы решить эту проблему, хороший подход должен был бы создать границу вокруг M делает его +2 размера в каждом измерении, и заселить, что с нулями, однако я действительно хотел бы, чтобы избежать изменения M, как мой код является более сложным, чем тот, который приведен в примере.

Я не могу найти способ сделать это, я немного застрял здесь. Любое предложение?

EDIT: @ Ответ работает, однако мне хотелось бы узнать, существует ли возможное решение с использованием этого метода линейной индексации.

% 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

Является ли этот линейный индексный материал существенным? Потому что это довольно легко обращаться с граничными условиями используется индекс индексации и min и max так:

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]) 

Таким образом, вы можете также легко обобщить это, если вы хотите более широкий район:

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]) 

или если вы любили держать форму куба, то возможно:

neighbours = neighbourhood; 
neighbours(mid) = NaN; 

Если вы хотите использовать это много раз в вашем коде это, вероятно, лучше реорганизовать его в качестве функции т-файла, который просто возвращает индексы:

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 

Answer 2

Базовая теория: Продлить входной массив left-right, up-down, one more on each sides of the third dimension с NaNs. Это позволило бы нам использовать единую сетку 3x3x3, а затем использовать эти NaNs для обнаружения элементов, выходящих за границы входного массива, и как таковые должны быть отброшены.

код

%// Initializations 
sz_ext = size(M)+2; %// Get size of padded/extended input 3D array 
M_ext = NaN(sz_ext); %// Initialize extended array 
M_ext(2:end-1,2:end-1,2:end-1) = M; %// Insert values from M into it 

%// Important stuff here : Calculate linear offset indices within one 3D slice 
%// then for neighboring 3D slices too 
offset2D = bsxfun(@plus,[-1:1]',[-1:1]*sz_ext(1)); %//' 
offset3D = bsxfun(@plus,offset2D,permute([-1:1]*sz_ext(1)*sz_ext(2),[1 3 2])); 

%// Get linear indices for all points 
points_linear_idx = sub2ind(size(M_ext),point(:,1)+1,point(:,2)+1,point(:,3)+1); 
%// Linear indices for all neighboring elements for all points; index into M_ext 
neigh26 = M_ext(bsxfun(@plus,offset3D,permute(points_linear_idx,[4 3 2 1]))) 

Как использовать: Таким образом, каждый срез в 4th измерения представляют собой 27 элементов (соседние плюс сам элемент) в виде 3x3x3 массива.Следовательно, neigh26 будет массивом 3x3x3xN, где N - это количество точек в массиве точек.

Пример: В качестве примера, давайте предположим, что некоторые случайные значения в M и Point -

M=rand(11,11,11); 
point = [ 
    1 1 4; 
    1 7 1] 

На запуск ранее кода с этими входами, я получаю что-то вроде этого -

neigh26(:,:,1,1) = 
     NaN  NaN  NaN 
     NaN 0.5859 0.4917 
     NaN 0.6733 0.6688 
neigh26(:,:,2,1) = 
     NaN  NaN  NaN 
     NaN 0.0663 0.5544 
     NaN 0.3440 0.3664 
neigh26(:,:,3,1) = 
     NaN  NaN  NaN 
     NaN 0.3555 0.1257 
     NaN 0.4424 0.9577 
neigh26(:,:,1,2) = 
    NaN NaN NaN 
    NaN NaN NaN 
    NaN NaN NaN 
neigh26(:,:,2,2) = 
     NaN  NaN  NaN 
    0.7708 0.3712 0.2866 
    0.7088 0.3743 0.2326 
neigh26(:,:,3,2) = 
     NaN  NaN  NaN 
    0.4938 0.5051 0.9416 
    0.1966 0.0213 0.8036