Milan_gvf.m
上传用户:zlding2008
上传日期:2013-05-13
资源大小:1914k
文件大小:5k
源码类别:

2D图形编程

开发平台:

Matlab

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