THESE/codes/snake/basic_code/Snake2D.m

   1 function [P,J]=Snake2D(I,P,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 % [O,J]=Snake2D(I,P,Options)\r
   8 %  \r
   9 % inputs,\r
  10 %   I : An Image of type double preferable ranged [0..1]\r
  11 %   P : List with coordinates descriping the rough contour N x 2\r
  12 %   Options : A struct with all snake options\r
  13 %   \r
  14 % outputs,\r
  15 %   O : List with coordinates of the final contour M x 2\r
  16 %   J : Binary image with the segmented region\r
  17 %\r
  18 % options (general),\r
  19 %  Option.Verbose : If true show important images, default false\r
  20 %  Options.nPoints : Number of contour points, default 100\r
  21 %  Options.Gamma : Time step, default 1\r
  22 %  Options.Iterations : Number of iterations, default 100\r
  23 %\r
  24 % options (Image Edge Energy / Image force))\r
  25 %  Options.Sigma1 : Sigma used to calculate image derivatives, default 10\r
  26 %  Options.Wline : Attraction to lines, if negative to black lines otherwise white\r
  27 %                    lines , default 0.04\r
  28 %  Options.Wedge : Attraction to edges, default 2.0\r
  29 %  Options.Wterm : Attraction to terminations of lines (end points) and\r
  30 %                    corners, default 0.01\r
  31 %  Options.Sigma2 : Sigma used to calculate the gradient of the edge energy\r
  32 %                    image (which gives the image force), default 20\r
  33 %\r
  34 % options (Gradient Vector Flow)\r
  35 %  Options.Mu : Trade of between real edge vectors, and noise vectors,\r
  36 %                default 0.2. (Warning setting this to high >0.5 gives\r
  37 %                an instable Vector Flow)\r
  38 %  Options.GIterations : Number of GVF iterations, default 0\r
  39 %  Options.Sigma3 : Sigma used to calculate the laplacian in GVF, default 1.0\r
  40 %\r
  41 % options (Snake)\r
  42 %  Options.Alpha : Membrame energy  (first order), default 0.2\r
  43 %  Options.Beta : Thin plate energy (second order), default 0.2\r
  44 %  Options.Delta : Baloon force, default 0.1\r
  45 %  Options.Kappa : Weight of external image force, default 2\r
  46 %\r
  47 %\r
  48 % Literature:\r
  49 %   - Michael Kass, Andrew Witkin and Demetri TerzoPoulos "Snakes : Active\r
  50 %       Contour Models", 1987\r
  51 %   - Jim Ivins amd John Porrill, "Everything you always wanted to know\r
  52 %       about snakes (but wer afraid to ask)\r
  53 %   - Chenyang Xu and Jerry L. Prince, "Gradient Vector Flow: A New\r
  54 %       external force for Snakes\r
  55 %\r
  56 % Example, Basic:\r
  57 %\r
  58 %  % Read an image\r
  59 %   I = imread('testimage.png');\r
  60 %  % Convert the image to double data type\r
  61 %   I = im2double(I); \r
  62 %  % Show the image and select some points with the mouse (at least 4)\r
  63 %  %figure, imshow(I); [y,x] = getpts;\r
  64 %   y=[182 233 251 205 169];\r
  65 %   x=[163 166 207 248 210];\r
  66 %  % Make an array with the clicked coordinates\r
  67 %   P=[x(:) y(:)];\r
  68 %  % Start Snake Process\r
  69 %   Options=struct;\r
  70 %   Options.Verbose=true;\r
  71 %   Options.Iterations=300;\r
  72 %   [O,J]=Snake2D(I,P,Options);\r
  73 %  % Show the result\r
  74 %   Irgb(:,:,1)=I;\r
  75 %   Irgb(:,:,2)=I;\r
  76 %   Irgb(:,:,3)=J;\r
  77 %   figure, imshow(Irgb,[]); \r
  78 %   hold on; plot([O(:,2);O(1,2)],[O(:,1);O(1,1)]);\r
  79 %  \r
  80 % Example, GVF:\r
  81 %   I=im2double(imread('testimage2.png'));\r
  82 %   x=[96 51 98 202 272 280 182];\r
  83 %   y=[63 147 242 262 211 97 59];\r
  84 %   P=[x(:) y(:)];\r
  85 %   Options=struct;\r
  86 %   Options.Verbose=true;\r
  87 %   Options.Iterations=400;\r
  88 %   Options.Wedge=2;\r
  89 %   Options.Wline=0;\r
  90 %   Options.Wterm=0;\r
  91 %   Options.Kappa=4;\r
  92 %   Options.Sigma1=8;\r
  93 %   Options.Sigma2=8;\r
  94 %   Options.Alpha=0.1;\r
  95 %   Options.Beta=0.1;\r
  96 %   Options.Mu=0.2;\r
  97 %   Options.Delta=-0.1;\r
  98 %   Options.GIterations=600;\r
  99 %   [O,J]=Snake2D(I,P,Options);\r
 100 %   \r
 101 % Function is written by D.Kroon University of Twente (July 2010)\r
 102 \r
 103 % Process inputs\r
 104 defaultoptions=struct('Verbose',false,'nPoints',100,'Wline',0.04,'Wedge',2,'Wterm',0.01,'Sigma1',10,'Sigma2',20,'Alpha',0.2,'Beta',0.2,'Delta',0.1,'Gamma',1,'Kappa',2,'Iterations',100,'GIterations',0,'Mu',0.2,'Sigma3',1);\r
 105 if(~exist('Options','var')), \r
 106     Options=defaultoptions; \r
 107 else\r
 108     tags = fieldnames(defaultoptions);\r
 109     for i=1:length(tags)\r
 110          if(~isfield(Options,tags{i})), Options.(tags{i})=defaultoptions.(tags{i}); end\r
 111     end\r
 112     if(length(tags)~=length(fieldnames(Options))), \r
 113         warning('snake:unknownoption','unknown options found');\r
 114     end\r
 115 end\r
 116 \r
 117 % Convert input to double\r
 118 I = double(I);\r
 119 \r
 120 % If color image convert to grayscale\r
 121 if(size(I,3)==3), I=rgb2gray(I); end\r
 122 \r
 123 % The contour must always be clockwise (because of the balloon force)\r
 124 P=MakeContourClockwise2D(P);\r
 125 \r
 126 % Make an uniform sampled contour description\r
 127 P=InterpolateContourPoints2D(P,Options.nPoints);\r
 128 \r
 129 % Transform the Image into an External Energy Image\r
 130 Eext = ExternalForceImage2D(I,Options.Wline, Options.Wedge, Options.Wterm,Options.Sigma1);\r
 131 \r
 132 % Make the external force (flow) field.\r
 133 Fx=ImageDerivatives2D(Eext,Options.Sigma2,'x');\r
 134 Fy=ImageDerivatives2D(Eext,Options.Sigma2,'y');\r
 135 Fext(:,:,1)=-Fx*2*Options.Sigma2^2;\r
 136 Fext(:,:,2)=-Fy*2*Options.Sigma2^2;\r
 137 \r
 138 % Do Gradient vector flow, optimalization\r
 139 Fext=GVFOptimizeImageForces2D(Fext, Options.Mu, Options.GIterations, Options.Sigma3);\r
 140 \r
 141 % Show the image, contour and force field\r
 142 if(Options.Verbose)\r
 143     h=figure; set(h,'render','opengl')\r
 144      subplot(2,2,1),\r
 145       imshow(I,[]); \r
 146       hold on; plot(P(:,2),P(:,1),'b.'); hold on;\r
 147       title('The image with initial contour')\r
 148      subplot(2,2,2),\r
 149       imshow(Eext,[]); \r
 150       title('The external energy');\r
 151      subplot(2,2,3), \r
 152       [x,y]=ndgrid(1:10:size(Fext,1),1:10:size(Fext,2));\r
 153       imshow(I), hold on; quiver(y,x,Fext(1:10:end,1:10:end,2),Fext(1:10:end,1:10:end,1));\r
 154       title('The external force field ')\r
 155      subplot(2,2,4), \r
 156       imshow(I), hold on; plot(P(:,2),P(:,1),'b.'); \r
 157       title('Snake movement ')\r
 158     drawnow\r
 159 end\r
 160 \r
 161 \r
 162 % Make the interal force matrix, which constrains the moving points to a\r
 163 % smooth contour\r
 164 S=SnakeInternalForceMatrix2D(Options.nPoints,Options.Alpha,Options.Beta,Options.Gamma);\r
 165 h=[];\r
 166 for i=1:Options.Iterations\r
 167     P=SnakeMoveIteration2D(S,P,Fext,Options.Gamma,Options.Kappa,Options.Delta);\r
 168 \r
 169     % Show current contour\r
 170     if(Options.Verbose)\r
 171         if(ishandle(h)), delete(h), end\r
 172         h=plot(P(:,2),P(:,1),'r.');\r
 173         c=i/Options.Iterations;\r
 174         plot([P(:,2);P(1,2)],[P(:,1);P(1,1)],'-','Color',[c 1-c 0]);  drawnow\r
 175     end\r
 176 end\r
 177 \r
 178 if(nargout>1)\r
 179      J=DrawSegmentedArea2D(P,size(I));\r
 180 end\r
 181 \r