MATLAB identify adjacient regions in 3D image

Question

I have a 3D image, divided into contiguous regions where each voxel has the same value. The value assigned to this region is unique to the region and serves as a label. The example image below describes the 2D case:

     1 1 1 1 2 2 2
     1 1 1 2 2 2 3
Im = 1 4 1 2 2 3 3
     4 4 4 4 3 3 3
     4 4 4 4 3 3 3

I want to create a graph describing adjaciency between these regions. In the above case, this would be:

    0 1 0 1
A = 1 0 1 1
    0 1 0 1
    1 1 1 0

I'm looking for a speedy solution to do this for large 3D images in MATLAB. I came up with a solution that iterates over all regions, which takes 0.05s per iteration - unfortunately, this will take over half an hour for an image with 32'000 regions. Does anybody now a more elegant way of doing this? I'm posting the current algorithm below:

labels = unique(Im); % assuming labels go continuously from 1 to N
A = zeros(labels);

for ii=labels
  % border mask to find neighbourhood
  dil = imdilate( Im==ii, ones(3,3,3) );
  border = dil - (Im==ii);

  neighLabels = unique( Im(border>0) );
  A(ii,neighLabels) = 1;
end

imdilate is the bottleneck I would like to avoid.

Thank you for your help!

Answer 1

I came up with a solution which is a combination of Divakar's and teng's answers, as well as my own modifications and I generalised it to the 2D or 3D case.

To make it more efficient, I should probably pre-allocate the r and c, but in the meantime, this is the runtime:

• For a 3D image of dimension 117x159x126 and 32000 separate regions: 0.79s
• For the above 2D example: 0.004671s with this solution, 0.002136s with Divakar's solution, 0.03995s with teng's solution.

I haven't tried extending the winner (Divakar) to the 3D case, though!

noDims = length(size(Im));
validim = ones(size(Im))>0;
labels = unique(Im);

if noDims == 3
    Im = padarray(Im,[1 1 1],'replicate', 'post');
    shifts = {[-1 0 0] [0 -1 0] [0 0 -1]};
elseif noDims == 2
    Im = padarray(Im,[1 1],'replicate', 'post');
    shifts = {[-1 0] [0 -1]};
end

% get value of the neighbors for each pixel
% by shifting the image in each direction
r=[]; c=[];
for i = 1:numel(shifts)
    tmp = circshift(Im,shifts{i});
    r = [r ; Im(validim)];
    c = [c ; tmp(validim)]; 
end

A = sparse(r,c,ones(size(r)), numel(labels), numel(labels) );
% make symmetric, delete diagonal
A = (A+A')>0;
A(1:size(A,1)+1:end)=0;

Thanks for the help!