THESE/codes/snake/basic_code/GVFOptimizeImageForces3D.m

   1 function Fext=GVFOptimizeImageForces3D(Fext, Mu, Iterations, Sigma)\r
   2 % This function "GVFOptimizeImageForces" does gradient vector flow (GVF)\r
   3 % on a vector field. GVF gives the edge-forces a larger capature range,\r
   4 % to make the snake also reach concave regions\r
   5 %\r
   6 % Fext = GVFOptimizeImageForces3D(Fext, Mu, Iterations, Sigma) \r
   7 % \r
   8 % inputs,\r
   9 %   Fext : The image force vector field N x M x 3\r
  10 %   Mu : Is a trade of scalar between noise and real edge forces\r
  11 %   Iterations : The number of GVF itterations\r
  12 %   Sigma : Used when calculating the Laplacian\r
  13 % \r
  14 % outputs,\r
  15 %   Fext : The GVF optimized image force vector field\r
  16 %\r
  17 % Function is written by D.Kroon University of Twente (July 2010)\r
  18 \r
  19 % Squared magnitude of force field\r
  20 Fx= Fext(:,:,:,1);\r
  21 Fy= Fext(:,:,:,2);\r
  22 Fz= Fext(:,:,:,3);\r
  23 \r
  24 % Calculate magnitude\r
  25 sMag = Fx.^2+ Fy.^2 + Fz.^2;\r
  26 \r
  27 % Set new vector-field to initial field\r
  28 u=Fx; v=Fy; w=Fz;\r
  29 \r
  30 % Iteratively perform the Gradient Vector Flow (GVF)\r
  31 for i=1:Iterations,\r
  32   % Calculate Laplacian\r
  33   Uxx=ImageDerivatives3D(u,Sigma,'xx');\r
  34   Uyy=ImageDerivatives3D(u,Sigma,'yy');\r
  35   Uzz=ImageDerivatives3D(u,Sigma,'zz');\r
  36 \r
  37   % Update the vector field\r
  38   u = u + Mu*(Uxx+Uyy+Uzz) - sMag.*(u-Fx);\r
  39   clear('Uxx','Uyy','Uzz');\r
  40   \r
  41   Vxx=ImageDerivatives3D(v,Sigma,'xx');\r
  42   Vyy=ImageDerivatives3D(v,Sigma,'yy');\r
  43   Vzz=ImageDerivatives3D(v,Sigma,'zz');\r
  44    \r
  45   v = v + Mu*(Vxx+Vyy+Vzz) - sMag.*(v-Fy);\r
  46   clear('Vxx','Vyy','Vzz');\r
  47   \r
  48   Wxx=ImageDerivatives3D(w,Sigma,'xx');\r
  49   Wyy=ImageDerivatives3D(w,Sigma,'yy');\r
  50   Wzz=ImageDerivatives3D(w,Sigma,'zz');\r
  51   \r
  52   w = w + Mu*(Wxx+Wyy+Wzz) - sMag.*(w-Fz);\r
  53   clear('Wxx','Wyy','Wzz');\r
  54 end\r
  55 \r
  56 Fext(:,:,:,1) = u;\r
  57 Fext(:,:,:,2) = v;\r
  58 Fext(:,:,:,3) = w;