MATLAB 对离散点进行圆拟合

function [z r residual] = fitcircle（x varargin）
%FITCIRCLE    least squares circle fit
%   
%   [Z R] = FITCIRCLE（X） fits a circle to the N points in X minimising
%   geometric error （sum of squared distances from the points to the fitted
%   circle） using nonlinear least squares （Gauss Newton）
%       Input
%           X : 2xN array of N 2D points with N >= 3
%       Output
%           Z : center of the fitted circle
%           R : radius of the fitted circle
%
%   [Z R] = FITCIRCLE（X ‘linear‘） fits a circle using linear least
%   squares minimising the algebraic error （residual from fitting system
%   of the form ax‘x + b‘x + c = 0）
%
%   [Z R] = FITCIRCLE（X Property Value ...） allows parameters to be
%   passed to the internal Gauss Newton method. Property names can be
%   supplied as any unambiguous contraction of the property name and are 
%   case insensitive e.g. FITCIRCLE（X ‘t‘ 1e-4） is equivalent to
%   FITCIRCLE（X ‘tol‘ 1e-4）. Valid properties are:
%
%       Property:                 Value:
%       --------------------------------
%       maxits                    positive integer default 100
%           Sets the maximum number of iterations of the Gauss Newton
%           method
%
%       tol                       positive constant default 1e-5
%           Gauss Newton converges when the relative change in the solution
%           is less than tol
%
%   [X R RES] = fitcircle（...） returns the 2 norm of the residual from 
%   the least squares fit
%
%   Example:
%       x = [1 2 5 7 9 3; 7 6 8 7 5 7];
%       % Get linear least squares fit
%       [zl rl] = fitcircle（x ‘linear‘）
%       % Get true best fit
%       [z r] = fitcircle（x）
%
%   Reference: “Least-squares fitting of circles and ellipses“ W. Gander
%   G. Golub R. Strebel - BIT Numerical Mathematics 1994 Springer

% This implementation copyright Richard Brown 2007 but is freely
% available to copy use or modify as long as this line is maintained

error（nargchk（1 5 nargin ‘struct‘））

% Default parameters for Gauss Newton minimisation
params.maxits = 100;
params.tol    = 1e-5;

% Check x and get user supplied parameters
[x fNonlinear params] = parseinputs（x params varargin{:}）;

% Convenience variables
m  = size（x 2）;
x1 = x（1 :）‘;
x2 = x（2 :）‘;


% 1） Compute best fit w.r.t. algebraic error using linear least squares
% 
% Circle is represented as a matrix quadratic form
%     ax‘x + b‘x + c = 0
% Linear least squares estimate found by minimising Bu = 0 s.t. norm（u） = 1
%     where u = [a; b; c]

% Form the coefficient matrix
B = [x1.^2 + x2.^2 x1 x2 ones（m 1）];

% Least squares estimate is right singular vector corresp. to smallest
% singular value of B
[U S V] = svd（B）;
u = V（: 4）;

% For clarity set the quadratic form variables
a = u（1）;
b = u（2:3）;
c = u（4）;

% Convert to centre/radius
z = -b / （2*a）;
r = sqrt（（norm（b）/（2*a））^2 - c/a）;

% 2