pare_varpro.m

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

function [z, a, b, alpha, phi, step] = ...
         pare_varpro (X, z, a, b, alpha, show);
%PARE_VARPRO    Geometric ellipse fit using varpro
%
% [z, a, b, alpha, phi, step] = ...
%      pare_varpro (X, z, a, b, alpha, show{0});
% computes the best fit ellipse in parameterform
% x = z(1) + a  cos(phi-alpha), y = z(2) + b sin(phi-alpha)
% using the varpro algorithm.
%
% X: given points <X(i,1), X(i,2)>
% z, a, b, alpha: starting values
% show: if (show), test output
%
% z, a, b, alpha: ellipse found
% phi: values for the nearest points (in parametric form)
% step: nof iterations

  m = size(X,1);

  s = sin(alpha);
  c = cos(alpha);
  Q = [c -s; s c];

  % compute initial approximations for phi_i
  du  = Q'*[X(:,1)-z(1) X(:,2)-z(2)]';
  phi = angle(du(1,:)/a + sqrt(-1)*du(2,:)/b)';

  [vn, vl, err, step] = varpro ('pare_varpro_hook', ...
                          [X(:,1);X(:,2)], [phi;alpha], [a;b;z]);

  phi   = vn(1:m);
  alpha = vn(1+m);
  a     = vl(1);
  b     = vl(2);
  z     = vl(3:4);

end % pare_varpro