algellipse.m

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

function [z, a, b, alpha, err] = algellipse (X, W, show);
%ALGELLIPSE     Algebraic least square ellipse fit
%
%       [z, a, b, alpha, err] = ...
%         algellipse (X, W{default ones}, show{default 0})
%       fits an ellipse by minimizing the "algebraic distance"
%       in the least squares sense x'A x + b'x + c = 0
%       weighting the i-th data by W(i)
%
%       X: given points Pi = [X(i,1), X(i,2)]
%       W: weight W(i) for the i-th equation
%       show: if (show == 1) make test output
%
%       z, a, b, alpha: parameters for found ellipse
%       err: error indication
%         if (err == 1), not an ellipse
%         if (err == 0), ok

  if (nargin < 2), W = ones(size(X,1), 1); end;
  if (nargin < 3), show = 0; end;

  [U S V] = svd(diag(W) * [X(:,1).^2  X(:,1).*X(:,2) X(:,2).^2 ...
                          X(:,1) X(:,2) ones(size(X(:,1)))]);
  u   = V(:,6);

  [z, a, b, alpha, err] = ellipse_params (u, show);

end % algellipse