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

2D图形编程

开发平台:

Matlab

  1. function [x,y] = snakedeform(x,y,alpha,beta,gamma,kappa, kappap, fx,fy,ITER)
  2. % SNAKEDEFORM2  Deform snake in the given external force field with
  3. %     pressure force added
  4. %     [x,y] = snakedeform(x,y,alpha,beta,gamma,kappa,kappap,fx,fy,ITER)
  5. %
  6. %     alpha:   elasticity parameter
  7. %     beta:    rigidity parameter
  8. %     gamma:   viscosity parameter
  9. %     kappa:   external force weight
  10. %     kappap:  pressure force weight
  11. %     fx,fy:   external force field
  12. %      Chenyang Xu and Jerry L. Prince, 4/1/95, 6/17/97
  13. %      Copyright (c) 1995-97 by Chenyang Xu and Jerry L. Prince
  14. %      Image Analysis and Communications Lab, Johns Hopkins University
  15. % generates the parameters for snake
  16. N = length(x);
  17. alpha = alpha* ones(1,N); 
  18. beta = beta*ones(1,N);
  19. % produce the five diagnal vectors
  20. alpham1 = [alpha(2:N) alpha(1)];
  21. alphap1 = [alpha(N) alpha(1:N-1)];
  22. betam1 = [beta(2:N) beta(1)];
  23. betap1 = [beta(N) beta(1:N-1)];
  24. a = betam1;
  25. b = -alpha - 2*beta - 2*betam1;
  26. c = alpha + alphap1 +betam1 + 4*beta + betap1;
  27. d = -alphap1 - 2*beta - 2*betap1;
  28. e = betap1;
  29. % generate the parameters matrix
  30. A = diag(a(1:N-2),-2) + diag(a(N-1:N),N-2);
  31. A = A + diag(b(1:N-1),-1) + diag(b(N), N-1);
  32. A = A + diag(c);
  33. A = A + diag(d(1:N-1),1) + diag(d(N),-(N-1));
  34. A = A + diag(e(1:N-2),2) + diag(e(N-1:N),-(N-2));
  35. invAI = inv(A + gamma * diag(ones(1,N)));
  36. for count = 1:ITER,
  37.    vfx = interp2(fx,x,y,'*linear');
  38.    vfy = interp2(fy,x,y,'*linear');
  39.   
  40.    % adding the pressure force 
  41.    xp = [x(2:N);x(1)];    yp = [y(2:N);y(1)]; 
  42.    xm = [x(N);x(1:N-1)];  ym = [y(N);y(1:N-1)]; 
  43.    
  44.    qx = xp-xm;  qy = yp-ym;
  45.    pmag = sqrt(qx.*qx+qy.*qy);
  46.    px = qy./pmag;     py = -qx./pmag;
  47.    % deform snake
  48.    x = invAI * (gamma* x + kappa*vfx + kappap.*px);
  49.    y = invAI * (gamma* y + kappa*vfy + kappap.*py);
  50. end