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!

Answer 2

Try this out -

Im = padarray(Im,[1 1],'replicate');

labels = unique(Im);
box1 = [-size(Im,1)-1 -size(Im,1) -size(Im,1)+1 -1 1 size(Im,1)-1 size(Im,1) size(Im,1)+1];

mat1 = NaN(numel(labels),numel(labels));
for k2=1:numel(labels)
    a1 = find(Im==k2);
    for k1=1:numel(labels)
        a2 = find(Im==k1);
        t1 = bsxfun(@plus,a1,box1);
        t2 = bsxfun(@eq,t1,permute(a2,[3 2 1]));
        mat1(k2,k1) = any(t2(:));
    end
end
mat1(1:size(mat1,1)+1:end)=0;

If it works for you, share with us the runtimes as comparison? Would love to see if the coffee brews any faster than half an hour!

Answer 3

Below is my attempt.

 Im = [1 1 1 1 2 2 2;
     1 1 1 2 2 2 3;
     1 4 1 2 2 3 3;
     4 4 4 4 3 3 3;
     4 4 4 4 3 3 3];

 % mark the borders
 validim = zeros(size(Im));
 validim(2:end-1,2:end-1) = 1;

 % get value of the 4-neighbors for each pixel
 % by shifting the images 4 times in each direction
 numNeighbors = 4;
 adj = zeros([prod(size(Im)),numNeighbors]);
 shifts = {[0 1] [0 -1] [1 0] [-1 0]};
 for i = 1:numNeighbors
     tmp = circshift(Im,shifts{i});
     tmp(validim == 0) = nan;
     adj(:,i) = tmp(:);
 end

 % mark neighbors where it does not eq Im
 imDuplicates = repmat(Im(:),[1 numNeighbors]);
 nonequals = adj ~= imDuplicates;
 % neglect the border
 nonequals(isnan(adj)) = 0;     
 % get these neighbor values and the corresponding Im value
 compared = [imDuplicates(nonequals == 1) adj(nonequals == 1)];

 % construct your 'A' % possibly could be more optimized here.
 labels = unique(Im);
 A = zeros(numel(labels));
 for i = 1:size(compared,1)
     A(compared(i,1),compared(i,2)) = 1;
 end

Answer 4

@Lisa Yours reasoning is elegant, though it obviously gives wrong answers for labels on the edges. Try this simple label matrix:

Im =

 1     2     2
 3     3     3
 3     4     4

The resulting adjacency matrix , according to your code is:

A =

 0     1     1     0
 1     0     1     1
 1     1     0     1
 0     1     1     0

which claims an adjacency between labels "2" and "4": obviously wrong. This happens simply because you are reading padded Im labels based on "validim" indices, which now doesn't match the new Im and goes all the way down to the lower borders.