Multi-dimensional version of bwboundaries

20 views (last 30 days)
Laurent
Laurent on 20 Aug 2013
Commented: Rik on 8 May 2020
Hello,
I am trying to get the boundaries of objects in a three-dimensional array. For 2D data, the bwboundaries function does this very well. Does a multi-dimensional version of this function exist?
I wrote something myself, by doing bwboundaries on each individual plane. The next steps are finding connected objects, renumbering them, calculating the new boundaries from this, and so on. This is all very slow on big arrays. So therefore I was hoping something, preferably faster, exists.
This is the code:
function [ B, C ] = bwboundaries2Dstack( BW )
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
B=cell(size(BW,3),1);
C=B;
N=zeros(size(BW,3),1);
for i=1:size(BW,3)
[B{i},C{i},N(i)]=bwboundaries(BW(:,:,i));
%first find holes and remove them
B{i}=B{i}(1:N(i));
for j=1:length(B{i})
B{i}{j}(:,3)=i;
end
C{i}(C{i}>N(i))=0;
end
%now make a 3D map of the objects found in 2D
maxC=0;
for i=2:size(BW,3)
connected=[];
maxC=max([maxC max(max(C{i-1}))]);
C{i}(C{i}>0)=C{i}(C{i}>0)+maxC; %make unique numbers compared to previous slice
tempim=zeros(size(BW,1),size(BW,2),2);
tempim(:,:,1)=C{i-1};
tempim(:,:,2)=C{i};
tempim=reshape(tempim,size(BW,1)*size(BW,2),2);
tempim=tempim(tempim(:,1)>0,:);
tempim=tempim(tempim(:,2)>0,:);
%for j=min(min(C{i}(C{i}>0))):max(max(C{i}))
%tempim=zeros(size(dataseg,1),size(dataseg,2));
%tempim(C{i}==j)=C{i-1}(C{i}==j);
if max(max(tempim))>0
overlaps=unique(tempim,'rows');
else
overlaps=[];
end
%now make a list showing which areas overlap from 2
%adjacent planes (plane i-1 and plane i)
if size(overlaps,2)>0
connected=[connected; overlaps];
end
%end
im1=C{i-1};
im2=C{i};
while ~isempty(connected)
currarea=connected(1,:);
keepind=[];
for j=2:size(connected,1) %find areas which are touching
if (sum(ismember(connected(j,1),currarea(:,1)))>0)||(sum(ismember(connected(j,2),currarea(:,2)))>0)
currarea=[currarea; connected(j,:)];
else
keepind=[keepind;j];
end
end
%now give every area the correct number
corrnumber=currarea(1,1); %all areas get the first number of the list
for j=2:size(currarea,1)
im1(C{i-1}==currarea(j,1))=corrnumber;
end
for j=1:size(currarea,1)
im2(C{i}==currarea(j,2))=corrnumber;
end
connected=connected(keepind,:);
end
C{i-1}=im1;
C{i}=im2;
%clear im1 im2
%i;
end
% convert C to a 3d matrix and renumber
Cnew=zeros(size(C{1},1),size(C{1},2),length(C));
for i=1:length(C)
Cnew(:,:,i)=C{i};
end
C=Cnew;
clear Cnew
allind=(sortrows(unique(C),1));
for i=2:length(allind)
C(C==allind(i))=i-1;
end
% Calculate the corresponding B
%clear B
B=cell(max(max(max(C))),1);
for i=1:max(max(max(C)))
%tempim=false(size(C));
%tempim(C==i)=1;
tempim=bsxfun(@eq,C,i);
%BW=bwperim(tempim);
zproj=sum(tempim,3);
xmin=find(sum(zproj,2)>0,1,'first');
xmax=find(sum(zproj,2)>0,1,'last');
ymin=find(sum(zproj,1)>0,1,'first');
ymax=find(sum(zproj,1)>0,1,'last');
planes=sum(sum(tempim,1),2);
planes=planes(:);
zmin=find(planes>0,1,'first');
zmax=find(planes>0,1,'last');
test=bwperim(tempim(xmin:xmax,ymin:ymax,zmin:zmax));
BW=false(size(tempim));
BW(xmin:xmax,ymin:ymax,zmin:zmax)=test;
[x,y,z]=ind2sub(size(BW),find(BW));
B{i}=[x y z];
end
end
The main problem are the for-loops, like this one:
for i=2:length(allind)
C(C==allind(i))=i-1;
end
and I don't see how I can speed this up.
  2 Comments
Matt Kindig
Matt Kindig on 20 Aug 2013
Edited: Matt Kindig on 20 Aug 2013
How would the boundary be defined for a 3D object? As a connected surface?
The bwconncomp() (on newer Matlab versions), or the older bwlabel() function, both work in 3D--does this help?
Laurent
Laurent on 20 Aug 2013
Edited: Laurent on 20 Aug 2013
Thank you very much. It seems that bwconncomp in combination with labelmatrix does the trick, I am trying it now. Not sure why I didn't find this myself...

Sign in to comment.

Answers (2)

Image Analyst
Image Analyst on 20 Aug 2013
Like Matt said, this is not so straightforward for a 3D surface. One thing you can do is to take your 3D binary object and erode it and subtract it from the original. That will get you a 3D volume where all the surface pixels are true and everything else is false
surfaceVoxels = binary3D - imerode(binary3D, true(3));
This will probably be good for your needs. If not, say what you want to do with it once you have it.
  5 Comments
Image Analyst
Image Analyst on 20 Aug 2013
Edited: Image Analyst on 20 Aug 2013
tempim = C>0; % Binarize.
But I really wonder why you need a list of all x,y,z voxels on the surface of each of your 3D blobs. Like I asked before, what do you want to do once you have that? It's possible that there is another way, like logical indexing.
Laurent
Laurent on 21 Aug 2013
Hi, thanks again for your comment and your suggestion. I have been looking into my code and I might be able to do what I want without getting this list of voxels.

Sign in to comment.


You Hao
You Hao on 8 May 2020
There is a much easier way:
d1=bwdist(bw1);
d2=bwdist(~bw1);
bwbd=(d1==d2+1);
  1 Comment
Rik
Rik on 8 May 2020
This doesn't provide the opportunity to use the 8 neighborhood instead of the 26 neighborhood. It also has two very expensive function calls, instead of a convolution (i.e. imerode).
How is this much easier?

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!