Understanding Voronoi Skeleton and extract the Algorithm ?

Question

mika 2014-2-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/114795-understanding-voronoi-skeleton-and-extract-the-algorithm

评论： mika 2014-2-18

Can any one help me to extract the Algorithm of the Voronoi Skeleton code below in order to understand it:

function [skel v e]=squelette(BW,varargin)
iptcheckinput(BW,{'numeric' 'logical'},{'real' 'nonsparse' '2d'}, ...
                mfilename, 'BW', 1);
  trim=0;
factor=-1;
boundary=false;
for i=1:length(varargin)
      switch lower(varargin{i})
          case 'trim'
              trim=varargin{i+1};
          case 'fast'
              factor=varargin{i+1};
          case 'boundary'
              boundary=true;
              if nargout~=2
                  error(['If boundary is specified, the function cannot 
  generate'...
                          char(10) 'the BW inage of the skeleton.' char(10)...
                          'Use [v e]=voronoiSkel(...) instead'])
              end
              if size(BW,2)~=2 
                  if size(BW,1)~=2
                      error('If you use the ''boundary'' option, the imput must 
be a 2 x n or n x 2 matrix')
                  else
                      BW=BW';
                  end
              end
      end
end
if trim<1;  trim=pi;  end
if factor<0; factor=1; end;
if ~islogical(BW) && ~boundary
      BW = (BW ~= 0);
end
%construct voronoi
if ~boundary
      b=bwboundaries(BW);
else
      b={BW};
end
if factor>1
      for i=1:length(b)
          inds=round(1:factor:size(b{i},1));
          b{i}=b{i}(inds,:);
      end
end
i=1;
inds=[];
while i<=length(b)
      if size(b{i},1)<4
          b(i)=[];
          continue;
      end
      inds(i)=length(b{i}); %#ok<AGROW>
      i=i+1;
end
inds=[0 cumsum(inds)];
p=cell2mat(b);
[v e]=costumVoronoi(p);
%clear bad vertices (bv) which are outside of the object.
if ~boundary
      rv=round(v);
      M=max(p);
      m=min(p);
      bv=v(:,1)<m(1)|v(:,1)>M(1)|v(:,2)>M(2)|v(:,2)<m(2);
      bv2=find(~bv);
      tmp=sub2ind(size(BW),rv(bv2,1),rv(bv2,2));
      bv2(BW(tmp))=[];
      bv=[find(bv); bv2];
else
      bv=find(~inpolygon(v(:,1),v(:,2),p(:,1),p(:,2)));
end
be=ismember(e(:,3),bv)|ismember(e(:,4),bv);
e=e(~be,:);
clear bv2 m M rv tmp;
% build distance table
D=cell(size(b));
for i=1:length(D);
      tmp=diff(b{i});
      tmp=sqrt(sum(tmp'.^2)');
      D{i}=[0 ;cumsum(tmp)];
end
% trim
be=false(size(e,1),1);
for i=1:size(e,1)
      i1=find(inds>=e(i,1),1,'first')-1;
      i2=find(inds>=e(i,2),1,'first')-1;
      if i1~=i2; continue; end;
      offset=inds(i1);
      contourDistance=abs(D{i1}(e(i,1)-offset)-D{i1}(e(i,2)-offset));
      contourDistance=min(contourDistance,D{i1}(end)-contourDistance);
      realDistance=norm(p(e(i,1),:)-p(e(i,2),:));
      if (contourDistance<realDistance*trim)
          be(i)=1;
      end
end
% keep only good edges
e=e(~be,3:4);
outputVertices=(nargout>1);
outputSkel=(nargout~=2);
if (outputSkel) 
      skel=false(size(BW));
      for i=1:size(e);
          v1=v(e(i,1),:);
          v2=v(e(i,2),:);
          t=linspace(0,1,max(ceil(1.3*norm(v2-v1)),4));
          x=v1(:,1).*t+(1-t).*v2(:,1);
          y=v1(:,2).*t+(1-t).*v2(:,2);
          inds=unique(round([x' y']),'rows');
          skel(sub2ind(size(skel),inds(:,1),inds(:,2)))=1;
      end
end
if outputVertices
      tmp=1:length(v);
      tmp=~ismember(tmp,e(:,1:2));
      inds=cumsum(tmp);
      e(:,1)=e(:,1)-inds(e(:,1))';
      e(:,2)=e(:,2)-inds(e(:,2))';
      v=v(~tmp,:);
end
if nargout==2
      skel=v;
      v=e;
end
end
function [v e]=costumVoronoi(V)
% Calculates the voronoi diagram of the vertices in V.
% v should be a n x 2 real matrix.
% qhull (www.qhull.org) MUST be executable from the current directory.
%
% Output: v is the matrix of the Voronoi vertices.
%         e is the matrix of the Voronoi edges. 
%            each row of e represents an edge in the following way: 
%            e(k,[1 2]) are the row indices (in V) of the points which
%            generated the k-th edge. 
%            e(k,[3 4]) are the row indices (in v) of the endpoints of the
%            k-th edge.
%         
% write to temp file
namein=tempname;
fid=fopen(namein,'w');
fprintf(fid,'%d\n%d\n',2,size(V,1));
fprintf(fid,'%d %d\n',V');
fclose(fid);
[a s]=dos(['qhull v p Fv TI ' namein]);
if (a~=0)
      delete(namein);
      error(['qhull returned an error:' char(10) s]);
end
[junk s]=strtok(s);
[Nv s]=strtok(s);
Nv=str2double(Nv);
[v pos]=textscan(s,'%f %f',Nv);
v=cell2mat(v);
s=s(pos:end);
[Ne s]=strtok(s);
e=double(cell2mat(textscan(s,'%d %d %d %d %d',Ne)));
e=e(:,2:5);
e(:,1:2)=e(:,1:2)+1;
e(e(:,3)==0,:)=[];
e(e(:,4)==0,:)=[];
%clean up
delete(namein);
end