Multi-dimensional version of bwboundaries
25 views (last 30 days)
Show older comments
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
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?
Answers (2)
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
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.
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
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?
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!