THESE/codes/snake/basic_code/Snake3D.m

   1 function FV=Snake3D(I,FV,Options)\r
   2 % This function SNAKE implements the basic snake segmentation. A snake is an \r
   3 % active (moving) contour, in which the points are attracted by edges and\r
   4 % other boundaries. To keep the contour smooth, an membrame and thin plate\r
   5 % energy is used as regularization.\r
   6 %\r
   7 % OV=Snake3D(I,FV,Options)\r
   8 %  \r
   9 % inputs,\r
  10 %   I : An Image of type double preferable ranged [0..1]\r
  11 %   FV : Structure with triangulated mesh, with list of faces FV.faces N x 3\r
  12 %        and list of vertices M x 3\r
  13 %   Options : A struct with all snake options\r
  14 %   \r
  15 % outputs,\r
  16 %   OV : Structure with triangulated mesh of the final surface\r
  17 %\r
  18 % options (general),\r
  19 %  Option.Verbose : If true show important images, default false\r
  20 %  Options.Gamma : Time step, default 1\r
  21 %  Options.Iterations : Number of iterations, default 100\r
  22 %\r
  23 % options (Image Edge Energy / Image force))\r
  24 %  Options.Sigma1 : Sigma used to calculate image derivatives, default 2\r
  25 %  Options.Wline : Attraction to lines, if negative to black lines otherwise white\r
  26 %                    lines , default 0.04\r
  27 %  Options.Wedge : Attraction to edges, default 2.0\r
  28 %  Options.Sigma2 : Sigma used to calculate the gradient of the edge energy\r
  29 %                    image (which gives the image force), default 2\r
  30 %\r
  31 % options (Gradient Vector Flow)\r
  32 %  Options.Mu : Trade of between real edge vectors, and noise vectors,\r
  33 %                default 0.2. (Warning setting this to high >0.5 gives\r
  34 %                an instable Vector Flow)\r
  35 %  Options.GIterations : Number of GVF iterations, default 0\r
  36 %  Options.Sigma3 : Sigma used to calculate the laplacian in GVF, default 1.0\r
  37 %\r
  38 % options (Snake)\r
  39 %  Options.Alpha : Membrame energy  (first order), default 0.2\r
  40 %  Options.Beta : Thin plate energy (second order), default 0.2\r
  41 %  Options.Delta : Baloon force, default 0.1\r
  42 %  Options.Kappa : Weight of external image force, default 2\r
  43 %  Options.Lambda : Weight which changes the direction of the image \r
  44 %                   potential force to the direction of the surface\r
  45 %                   normal, value default 0.8 (range [0..1])\r
  46 %                   (Keeps the surface from self intersecting)\r
  47 %\r
  48 % Literature:\r
  49 %   - Michael Kass, Andrew Witkin and Demetri TerzoPoulos "Snakes : Active\r
  50 %       Contour Models", 1987\r
  51 %   - Christoph Lurig, Leif Kobbelt, Thomas Ertl, "Hierachical solutions\r
  52 %       for the Deformable Surface Problem in Visualization"\r
  53 %\r
  54 % Example,\r
  55 %\r
  56 % load testvolume\r
  57 % load SphereMesh\r
  58 % Options=struct;\r
  59 % Options.Verbose=1;\r
  60 % Options.Wedge=0;\r
  61 % Options.Wline=-1;\r
  62 % Options.Alpha=0.2;\r
  63 % Options.Beta=0.2;\r
  64 % Options.Kappa=0.5;\r
  65 % Options.Delta=0.1000;\r
  66 % Options.Gamma=0.1000;\r
  67 % Options.Iterations=20;\r
  68 % Options.Sigma1=2;\r
  69 % Options.Sigma2=2;\r
  70 % Options.Lambda=0.8;\r
  71 % FV.vertices(:,1)=FV.vertices(:,1)+35;\r
  72 % FV.vertices(:,2)=FV.vertices(:,2)+25;\r
  73 % FV.vertices(:,3)=FV.vertices(:,3)+20;\r
  74 % OV=Snake3D(I,FV,Options)\r
  75 %\r
  76 % Function is written by D.Kroon University of Twente (July 2010)\r
  77 \r
  78 % Process inputs\r
  79 defaultoptions=struct('Verbose',false,'Wline',0.04,'Wedge',2,'Sigma1',2,'Sigma2',2,'Alpha',0.2,'Beta',0.2,'Delta',0.1,'Gamma',1,'Kappa',2,'Iterations',100,'GIterations',0,'Mu',0.2,'Sigma3',1,'Lambda',0.8);\r
  80 if(~exist('Options','var')), \r
  81     Options=defaultoptions; \r
  82 else\r
  83     tags = fieldnames(defaultoptions);\r
  84     for i=1:length(tags)\r
  85          if(~isfield(Options,tags{i})), Options.(tags{i})=defaultoptions.(tags{i}); end\r
  86     end\r
  87     if(length(tags)~=length(fieldnames(Options))), \r
  88         warning('snake:unknownoption','unknown options found');\r
  89     end\r
  90 end\r
  91 \r
  92 % Convert input to single if xintxx\r
  93 if(~strcmpi(class(I),'single')&&~strcmpi(class(I),'double'))\r
  94     I=single(I);\r
  95 end\r
  96 \r
  97 % The surface faces must always be clockwise (because of the balloon force)\r
  98 FV=MakeContourClockwise3D(FV);\r
  99 \r
 100 % Transform the Image into an External Energy Image\r
 101 Eext = ExternalForceImage3D(I,Options.Wline, Options.Wedge,Options.Sigma1);\r
 102 \r
 103 % Make the external force (flow) field.\r
 104 Fx=ImageDerivatives3D(Eext,Options.Sigma2,'x');\r
 105 Fy=ImageDerivatives3D(Eext,Options.Sigma2,'y');\r
 106 Fz=ImageDerivatives3D(Eext,Options.Sigma2,'z');\r
 107 \r
 108 Fext(:,:,:,1)=-Fx*2*Options.Sigma2^2;\r
 109 Fext(:,:,:,2)=-Fy*2*Options.Sigma2^2;\r
 110 Fext(:,:,:,3)=-Fz*2*Options.Sigma2^2;\r
 111 \r
 112 % Do Gradient vector flow, optimalization\r
 113 Fext=GVFOptimizeImageForces3D(Fext, Options.Mu, Options.GIterations, Options.Sigma3);\r
 114 \r
 115 % Show the image, contour and force field\r
 116 if(Options.Verbose)\r
 117      drawnow; pause(0.1);\r
 118      h=figure; set(h,'render','opengl')\r
 119      subplot(2,3,1),imshow(squeeze(Eext(:,:,round(end/2))),[]);\r
 120      subplot(2,3,2),imshow(squeeze(Eext(:,round(end/2),:)),[]);\r
 121      subplot(2,3,3),imshow(squeeze(Eext(round(end/2),:,:)),[]);\r
 122      subplot(2,3,4),imshow(squeeze(Fext(:,:,round(end/2),:))+0.5);\r
 123      subplot(2,3,5),imshow(squeeze(Fext(:,round(end/2),:,:))+0.5);\r
 124      subplot(2,3,6),imshow(squeeze(Fext(round(end/2),:,:,:))+0.5);\r
 125      h=figure; set(h,'render','opengl'); hold on;\r
 126      %patch(i,'facecolor',[1 0 0],'facealpha',0.1);\r
 127      ind=find(I(:)>0);\r
 128      [ix,iy,iz]=ind2sub(size(Eext),ind);\r
 129      plot3(ix,iy,iz,'b.');\r
 130      hold on;\r
 131      h=patch(FV,'facecolor',[1 0 0],'facealpha',0.1);\r
 132      drawnow; pause(0.1);\r
 133 end\r
 134 \r
 135 % Make the interal force matrix, which constrains the moving points to a\r
 136 % smooth contour\r
 137 S=SnakeInternalForceMatrix3D(FV,Options.Alpha,Options.Beta,Options.Gamma);\r
 138 for i=1:Options.Iterations\r
 139     FV=SnakeMoveIteration3D(S,FV,Fext,Options.Gamma,Options.Kappa,Options.Delta,Options.Lambda);\r
 140 \r
 141     % Show current contour\r
 142     if(Options.Verbose)\r
 143         if(ishandle(h));\r
 144             delete(h);\r
 145             h=patch(FV,'facecolor',[1 0 0],'facealpha',0.1);\r
 146             drawnow; %pause(0.1);\r
 147         end\r
 148     end\r
 149 end\r
 150 \r
 151 \r