THESE/codes/snake/gvf_dist_v4.2c/examples/gvf_ex.m

   1 % EXAMPLE     GVF snake examples on two simulated object boundaries.
   2 %
   3 % NOTE:
   4 %
   5 % traditional snake and distance snake differed from GVF snake only
   6 %   by using different external force field. In order to produce the
   7 %   corresponding external force field, use the following (all
   8 %   assuming edge map f is in memory).
   9 %
  10 % traditional snake:
  11 %   f0 = gaussianBlur(f,1);
  12 %   [px,py] = gradient(f0);
  13 %
  14 % distance snake:
  15 %   D = dt(f>0.5);  % distance transform (dt) require binary edge map
  16 %   [px,py] = gradient(-D);
  17 %
  18 % [px,py] is the external force field in both cases
  19 %
  20 % balloon model using a different matlab function "snakedeform2"
  21 % instead of "snakedeform". The external force could be any force
  22 % field generated above.
  23 %
  24 % an example of using it is as follows
  25 %       [x,y] = snakedeform2(x,y, alpha, beta, gamma, kappa, kappap, px,py,2);
  26 % do a "help snakedeform2" for more details
  27 %
  28
  29 %   Chenyang Xu and Jerry Prince 6/17/97
  30 %   Copyright (c) 1996-97 by Chenyang Xu and Jerry Prince
  31
  32    cd ..;   s = cd;   s = [s, '/snake']; path(s, path); cd examples;
  33
  34    help gvf_ex;
  35
  36    % ==== Example 1: U-shape object ====
  37
  38    % Read in the 64x64 U-shape image
  39      %[I,map] = rawread('../images/cochon64.pgm');
  40      I = imread('../images/cochon64.pgm');
  41    % Compute its edge map
  42      disp(' Compute edge map ...');
  43      f = gaussianBlur(1 - I/255,1);
  44
  45
  46    % Compute the GVF of the edge map f
  47      disp(' Compute GVF ...');
  48      [u,v] = GVF(f, 0.2, 80);
  49      disp(' Nomalizing the GVF external force ...');
  50      mag = sqrt(u.*u+v.*v);
  51      px = u./(mag+1e-10); py = v./(mag+1e-10);
  52
  53    % display the results
  54      figure(1);
  55      subplot(221); imdisp(I); title('test image');
  56      subplot(222); imdisp(f); title('edge map');
  57
  58      % display the gradient of the edge map
  59      [fx,fy] = gradient(f)
  60      subplot(223); quiver(fx,fy);
  61      axis off; axis equal; axis 'ij';     % fix the axis
  62      title('edge map gradient');
  63
  64      % display the GVF
  65      subplot(224); quiver(px,py);
  66      axis off; axis equal; axis 'ij';     % fix the axis
  67      title('normalized GVF field');
  68
  69    % snake deformation
  70      disp(' Press any key to start GVF snake deformation');
  71      pause;
  72      figure(1); subplot(221); cla;
  73      colormap(gray(64)); image(((1-f)+1)*40); axis('equal', 'off');
  74      t = 0:0.05:6.28;
  75      x = 32 + 30*cos(t);
  76      y = 32 + 30*sin(t);
  77      [x,y] = snakeinterp(x,y,2,0.5);
  78      snakedisp(x,y,'r')
  79      pause(1);
  80
  81      for i=1:50,
  82        [x,y] = snakedeform(x,y,2,0,1,0.1,px,py,5);
  83        [x,y] = snakeinterp(x,y,2,0.5);
  84        snakedisp(x,y,'r')
  85        title(['Deformation in progress,  iter = ' num2str(i*5)])
  86        pause(0.5);
  87      end
  88
  89      disp(' Press any key to display the final result');
  90      pause;
  91      cla;
  92      colormap(gray(64)); image(((1-f)+1)*40); axis('equal', 'off');
  93      snakedisp(x,y,'r')
  94      title(['Final result,  iter = ' num2str(i*5)]);
  95      disp(' ');
  96      disp(' Press any key to run the next example');
  97      pause;
  98
  99
 100    % ==== Example 2: Room like object ====
 101
 102    % Read in the 64x64 room image
 103      [I,map] = rawread('../images/room.pgm');
 104
 105    % Compute its edge map
 106      disp(' Compute edge map ...');
 107      f = 1 - I/255;
 108
 109    % Compute the GVF of the edge map f
 110      disp(' Compute GVF ...');
 111      [px,py] = GVF(f, 0.2, 80);
 112
 113    % display the results
 114      figure(1);
 115      subplot(221); imdisp(I); title('test image');
 116      subplot(222); imdisp(f); title('edge map');
 117
 118      % display the gradient of the edge map
 119      [fx,fy] = gradient(f);
 120      subplot(223); quiver(fx,fy);
 121      axis off; axis equal; axis 'ij';     % fix the axis
 122      title('edge map gradient');
 123
 124      % display the GVF
 125      subplot(224); quiver(px,py);
 126      axis off; axis equal; axis 'ij';     % fix the axis
 127      title('GVF field');
 128
 129    % snake deformation
 130      disp(' Press any key to start GVF snake deformation');
 131      pause;
 132      figure(1); subplot(221); cla;
 133      colormap(gray(64)); image(((1-f)+1)*40); axis('equal', 'off');
 134      t = 0:0.5:6.28;
 135      x = 32 + 3*cos(t);
 136      y = 32 + 3*sin(t);
 137      [x,y] = snakeinterp(x,y,2,0.5);
 138      snakedisp(x,y,'r');
 139      pause(1);
 140
 141      for i=1:15,
 142        [x,y] = snakedeform(x,y,0.05,0,1,2,px,py,5);
 143        [x,y] = snakeinterp(x,y,2,0.5);
 144        snakedisp(x,y,'r')
 145        title(['Deformation in progress,  iter = ' num2str(i*5)])
 146        pause(0.5);
 147      end
 148
 149      disp(' Press any key to display the final result');
 150      pause;
 151      figure(1); subplot(221); cla;
 152      colormap(gray(64)); image(((1-f)+1)*40); axis('equal', 'off');
 153      snakedisp(x,y,'r');
 154      title(['Final result,  iter = ' num2str(i*5)]);
 155
 156
 157
 158
 159
 160