ellipse_phi.m

来自「椭圆拟合的相关介绍与数学运算方法」· M 代码 · 共 47 行

function phi = ellipse_phi (X, a, b, phi, myeps);
%ELLIPSE_PHI    Compute nearest points to ellipse
%
% phi = ellipse_phi (X, a, b, phi{}, myeps{sqrt(myeps)});
% compute angles for nearest points on ellipse to given points X
% 
% X: given points <X(i,1),X(i,2)>
% a, b: ellipse parameters with a-axis in x-coordinate
% phi: starting values for iteration
% myeps: convergence limit
%
% phi: phi(i) angle of point on ellipse (in parametric form)
%      nearest to point <X(i,1), X(i,2)>

  if (nargin < 5), myeps = sqrt(eps); end;
  if (nargin < 4), phi = []; end;
  if (phi == []),  phi = angle (X(:,1)/a + sqrt(-1)*X(:,2)/b); end;
  
  m = size(X,1);
  for i = 1:m,
    par = [b^2 - a^2; a*X(i,1); b*X(i,2)];
    ni  = phi(i);
    step = 0;
    while (1),
      c = cos(ni);
      s = sin(ni);
      dni = - (par(1)*s*c + par(2)*s - par(3)*c) / ...
              (par(1)*(c^2 - s^2) + par(2)*c + par(3)*s);
      while ((X(i,1) - a*c)^2 + (X(i,2) - b*s)^2) < ...
            ((X(i,1) - a*cos(ni+dni))^2 + (X(i,2) - b*sin(ni+dni))^2),
        dni = -0.2*dni;
      end
      ni = ni + dni;
      if (abs(dni) <= sqrt(eps)) break; end;
      step = step + 1;
      if (step > 40), 
        disp ('warning: no convergence');
        str = sprintf ('angle, diff-angle = %g, %g', ni, dni);
        disp (str);
        break; 
      end
    end
    phi(i) = ni;
  end % for all points

end % ellipse_phi