Multi-dimensional version of bwboundaries

Question

Laurent 2013-8-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/85180-multi-dimensional-version-of-bwboundaries

评论： Rik 2020-5-8

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 个评论
显示无隐藏无

Matt Kindig 2013-8-20

编辑：Matt Kindig 2013-8-20

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 2013-8-20

编辑：Laurent 2013-8-20

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 2013-8-20

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/85180-multi-dimensional-version-of-bwboundaries#answer_94675

在 MATLAB Online 中打开

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 个评论
显示 3更早的评论隐藏 3更早的评论

Laurent 2013-8-20

在 MATLAB Online 中打开

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 2013-8-20

编辑：Image Analyst 2013-8-20

在 MATLAB Online 中打开

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 2013-8-21

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 2020-5-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/85180-multi-dimensional-version-of-bwboundaries#answer_431122

There is a much easier way:

d1=bwdist(bw1);

d2=bwdist(~bw1);

bwbd=(d1==d2+1);

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Rik 2020-5-8

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?

请先登录，再进行评论。

Multi-dimensional version of bwboundaries

2 个评论
显示无隐藏无

回答（2 个）

5 个评论
显示 3更早的评论隐藏 3更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

Multi-dimensional version of bwboundaries

2 个评论 显示 无隐藏 无

回答（2 个）

5 个评论 显示 3更早的评论隐藏 3更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

2 个评论
显示无隐藏无

5 个评论
显示 3更早的评论隐藏 3更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论