Multi-dimensional version of bwboundaries

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.

[Matt Kindig]

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]

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...

Answer 1

[Image Analyst]

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.

[Laurent]

Isn't that exactly what bwperim does?

[Image Analyst]

Yeah, you're right. I didn't know bwperim() worked on 3D images, but it does. So, will that work for you?

[Laurent]

Well, more or less. Thanks to Matt I have simplified and accelerated my code:

function [ B, C ] = bwboundaries2Dstackfast( BW )

CC=bwconncomp(BW);

C=labelmatrix(CC);

B=cell(max(max(max(C))),1);

% now find all the perimeters of each component in C.

for i=1:max(max(max(C)))

      tempim=bsxfun(@eq,C,i); %make binary image where the pixels of object i are true

      %make a ROI to speed up this process
      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');

      % now find the perimeter of object i
      test=bwperim(tempim(xmin:xmax,ymin:ymax,zmin:zmax));
      BW=false(size(tempim));
      BW(xmin:xmax,ymin:ymax,zmin:zmax)=test;

      % and create the coordinates of the border
      [x,y,z]=ind2sub(size(BW),find(BW));
      B{i}=[x y z];
  end

However the part where I determine the border/perimeter of each object is very slow, because my data contains many objects and my array is quite big (1000x1000x250 for example). Is there a way to speed it up?

According to the profiler the following lines are taking the most time:

tempim=bsxfun(@eq,C,i); 

zproj=sum(tempim,3);

[x,y,z]=ind2sub(size(BW),find(BW));

[Image Analyst]

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]

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.

Answer 2

[You Hao]

There is a much easier way:

d1=bwdist(bw1);

d2=bwdist(~bw1);

bwbd=(d1==d2+1);

[Rik]

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?