ransac fitline and fitellipse 实现Ransac拟合直线和椭圆

实现Ransac拟合直线和椭圆，可直接用。

function [apar1par2par3par4par5par6] = fitellipse（XY）

% FITELLIPSE Least-squares fit of ellipse to 2D points.
%        A = FITELLIPSE（XY） returns the parameters of the best-fit
%        ellipse to 2D points （XY）.
%        The returned vector A contains the center radii and orientation
%        of the ellipse stored as （Cx Cy Rx Ry theta_radians）
%
% Authors: Andrew Fitzgibbon Maurizio Pilu Bob Fisher
% Reference: “Direct Least Squares Fitting of Ellipses“ IEEE T-PAMI 1999
%
% This is a more bulletproof version than that in the paper incorporating
% scaling to reduce roundoff error correction of behaviour when the input
% data are on a perfect hyperbola and returns the geometric parameters
% of the ellipse rather than the coefficients of the quadratic form.
%
% Example: Run fitellipse without any arguments to get a demo
if nargin == 0
% Create an ellipse
t = linspace（02）;

Rx = 300
Ry = 200
Cx = 250
Cy = 150
Rotation = .4 % Radians

x = Rx * cos（t）;
y = Ry * sin（t）;
nx = x*cos（Rotation）-y*sin（Rotation） + Cx;
ny = x*sin（Rotation）+y*cos（Rotation） + Cy;
% Draw it
plot（nxny‘o‘）;
% Fit it
fitellipse（nxny）
% Note it returns （Rotation - pi/2） and swapped radii this is fine.
return
end

% normalize data
mx = mean（X）;
my = mean（Y）;
sx = （max（X）-min（X））/2;
sy = （max（Y）-min（Y））/2;

x = （X-mx）/sx;
y = （Y-my）/sy;

% Force to column vectors
x = x（:）;
y = y（:）;

% Build design matrix
D = [ x.*x x.*y y.*y x y ones（size（x）） ];

% Build scatter matrix
S = D‘*D;

% Build 6x6 constraint matrix
C（66） = 0; C（13） = -2; C（22） = 1; C（31） = -2;

% Solve eigensystem
[gevec geval] = eig（SC）;

% Find the negative eigenvalue
% [NegRNegC]=find（geval<0&~isinf（geval））;
% A=gevec（:NegC）;
I = find（real（diag（geval）） < 1e-8 & ~isinf（diag（geval）））;

% Extract eigenvector corresponding to negative eigenvalue
A = real（gevec（:I））;
%  m=（A‘*S*A）;
%  A=A*sqrt（geval（I）/m）;

% unnormalize
par = [
A（1）*sy*sy   ...
      A（2）*sx*sy   ...
      A（3）*sx*sx   ...
      -2*A（1）*sy*sy*mx - A（2）*sx*sy*my + A（4）*sx*sy*sy   ...
      -A（2）*sx*sy*mx - 2*A（3）*sx*sx*my + A（5）*sx*sx*sy   ...
      A（1）*sy*sy*mx*mx + A（2）*sx*sy*mx*my + A（3）*sx*sx*my*my   ...
      - A（4）*sx*sy*sy*mx - A（5）*sx*sx*sy*my   ...
      + A（6）*sx*sx*sy*sy   ...
      ]‘;

% Convert to geometric radii and centers

thetarad = 0.5*atan2（par（2）par（1） - par（3））;
cost = cos（thetarad）;
sint = sin（thetarad）;
sin_squared = sint.*sint;
cos_squared = cost.*cost;
cos_sin = sint .* cost;

Ao = par（6）;
Au =   par（4） .* cost + par（5） .* sint;
Av = - par（4） .* sint + par（5） .* cost;
Auu = par（1） .* cos_squared + par（3） .* sin_squared + par（2） .* cos_sin;
Avv = par（1） .* sin_squared + par（3） .* cos_squared - par（2） .* cos_sin;

% ROTATED = [Ao Au Av Auu Avv]

tuCentre = - Au./（2.*Auu）;
tvCentre = - Av./（2.*Avv）;
wCentre = Ao - Auu.*tuCentre.*tuCentre - Avv.*tvCentre.*tvCentre;

uCentre = tuCentre

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2012-06-21 14:25  ransac拟合代码\
     目录           0  2012-08-17 14:33  ransac拟合代码\椭圆拟合程序\
     文件        3553  2010-11-26 18:29  ransac拟合代码\椭圆拟合程序\fitellipse.m
     文件        4958  2010-11-26 19:35  ransac拟合代码\椭圆拟合程序\fitellipse1.m
     文件        1904  2010-11-26 18:34  ransac拟合代码\椭圆拟合程序\fitellipse2.m
     目录           0  2012-08-27 14:03  ransac拟合代码\直线拟合程序\
     文件     5019478  2011-12-28 17:07  ransac拟合代码\直线拟合程序\image4.bmp
     文件        1408  2012-05-04 16:44  ransac拟合代码\直线拟合程序\ransac.asv
     文件        1426  2012-08-09 16:00  ransac拟合代码\直线拟合程序\ransac.m
     文件      302586  2012-05-04 16:40  ransac拟合代码\直线拟合程序\硅臂直线.bmp
     文件      852282  2012-05-04 17:13  ransac拟合代码\直线拟合程序\硅臂直线2.bmp