0

I am a bit confused and would greatly appreciate some help.

I have read many posts about finding neighboring pixels, with this being extremely helpful:

http://blogs.mathworks.com/steve/2008/02/25/neighbor-indexing-2/

However I have trouble applying it on a 4D matrix (A) with size(A)=[8 340 340 15]. It represents 8 groups of 3D images (15 slices each) of which I want to get the neighbors. I am not sure which size to use in order to calculate the offsets. This is the code I tried, but I think it is not working because the offsets should be adapted for 4 dimensions? How can I do it without a loop?

%A is a 4D matrix with 0 or 1 values
Aidx = find(A); 

% loop here? 
[~,M,~,~] =size(A);
neighbor_offsets = [-1, M, 1, -M]';

neighbors_idx = bsxfun(@plus, Aidx', neighbor_offsets(:));
neighbors = B(neighbors_idx);

Thanks, ziggy

ziggy
  • 17
  • 4
  • Would something like bwlabeln work for this case? – ziggy Jun 23 '14 at 08:04
  • you have a bug: what happens if one of `Aidx` is at the image boundary? in that case your `neighbors` will point outside the image... – Shai Jun 23 '14 at 08:06
  • 1
    What kind of neighbors are you looking for? 4 neighbors along the second dimension? – Shai Jun 23 '14 at 08:06
  • I know for a fact in the specific datasets that I will not have Aidx in the image boundaries. I am looking for 4 neighbors on every slice on every group, so effectively on every (340*340) image. I will then need to make a mask of the nonzero neighbors and apply it on another matrix. – ziggy Jun 23 '14 at 08:12

3 Answers3

2

Not sure if I've understood your question but what about this sort of approach:

if you matrix is 1D:

M = rand(10,1);
N = M(k-1:k+1); %//immediate neighbours of k

However this could error if k is at the boundary. This is easy to fix using max and min:

N = M(max(k-1,1):min(k+1,size(M,1))

Now lets add a dimenion:

M = rand(10,10);
N = M(max(k1-1,1):min(k1+1,size(M,1), max(k2-1,1):min(k2+1,size(M,2))

That was easy, all you had to do was repeat the same index making the minor change of using size(M,2) for the boundary (and also I changed k to k1 and k2, you might find using an array for k instead of separate k1 and k2 variables works better i.e. k(1) and k(2))

OK so now lets skip to 4 dimensions:

M = rand(10,10,10,10);
N = M(max(k(1)-1,1):min(k(1)+1,size(M,1)), ...
      max(k(2)-1,1):min(k(2)+1,size(M,2)), ...
      max(k(3)-1,1):min(k(3)+1,size(M,3)), ...
      max(k(4)-1,1):min(k(4)+1,size(M,4)));  %// Also you can replace all the `size(M,i)` with `end` if you like

I know you said you didn't want a loop, but what about a really short loop just to refactor a bit and also make it generalized:

n=ndims(M);
ind{n} = 0;
for dim = 1:n
    ind{dim} = max(k(dim)-1,1):min(k(dim)+1,size(M,dim));
end
N = M(ind{:});
Dan
  • 45,079
  • 17
  • 88
  • 157
  • 1
    Thanks for this very detailed answer, Dan! Sorry, I don't have enough reputation to upvote it. I accepted Shai's answer because it is more compact and it works well. – ziggy Jun 23 '14 at 09:40
2

Have you considered using convn? 

msk = [0 1 0; 1 0 1; 0 1 0];
msk4d = permute( msk, [3 1 2 4] ); % make it 1-3-3-1 mask
neighbors_idx = find( convn( A, msk4d, 'same' ) > 0 );

You might find conndef useful for defining the basic msk in a general way.

Shai
  • 111,146
  • 38
  • 238
  • 371
  • What is `A` in this example? Is it the linear index of the element you're searching for the neighbours of? Surely you'd still need to give `convn` some information on the size of the original matrix? – Dan Jun 23 '14 at 08:29
  • @Dan `A` is the 4-D logical matrix as defined by the OP. – Shai Jun 23 '14 at 08:55
  • I like this solution too. Could you post the output to show that the neighbors are identified as well, please? – kkuilla Jun 23 '14 at 09:12
  • @kkuilla since I do not have the inputs - I cannot post the outputs ;) – Shai Jun 23 '14 at 12:43
  • Shai, one thing..in your second line with the permute function, the msk4d is not 1-3-3-1 mask, I think it is 1-3-3 – ziggy Jun 23 '14 at 18:44
  • @ziggy technically you are right, trailing singleton dimensions are not displayed. – Shai Jun 24 '14 at 06:11
0

Here's how to get the neighbors along the second dimension

sz = size( A );
ndims = numel(sz); % number of dimensions
[d{1:ndims}] = ind2sub( sz, find( A ) );
alongD = 2; % work along this dim
np = d{alongD} + 1;
sel = np <= sz( alongD ); % discard neighbors that fall outside image boundary
nm = d{alongD} - 1;
sel = sel & nm > 0; % discard neighbors that fall outside image boundary
d = cellfun( @(x) x(sel), d, 'uni', 0 ); 
neighbors = cat( 1, ...
                 ind2sub( sz, d{1:alongD-1}, np(sel), d{alongD+1:end} ),...
                 ind2sub( sz, d{1:alongD-1}, nm(sel), d{alongD+1:end} ) );
Shai
  • 111,146
  • 38
  • 238
  • 371