Home » Source Code » fisheye corrction » fitellipse.m

fitellipse.m ( File view )

  • By zhengzou 2014-05-09
  • View(s):31
  • Download(s):12
  • Point(s): 2
			function [Ru,uCentre vCentre] = fitellipse(X,Y)

% FITELLIPSE  Least-squares fit of ellipse to 2D points.
%        A = FITELLIPSE(X,Y) returns the parameters of the best-fit
%        ellipse to 2D points (X,Y).
%        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(0,2);
% 
%   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(nx,ny,'o');
%   % Fit it
%   fitellipse(nx,ny)
%   % 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(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;

% Solve eigensystem
[gevec, geval] = eig(S,C);

% Find the negative eigenvalue
I = find(real(diag(geval)) < 1e-8 & ~isinf(diag(geval)));

% Extract eigenvector corresponding to negative eigenvalue
A = real(gevec(:,I));

% 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 .* cost - tvCentre .* sint;
vCentre = tuCentre .* sint + tvCentre .* cost;

Ru = -wCentre./Auu;
Rv = -wCentre./Avv;

Ru = sqrt(abs(Ru)).*sign(Ru);
Rv = sqrt(abs(Rv)).*sign(Rv);

%a = [uCentre, vCentre, Ru, Rv, thetarad];

			
...
Expand> <Close

Want complete source code? Download it here

Point(s): 2

Download
0 lines left, continue to read
Sponsored links

File list

Tips: You can preview the content of files by clicking file names^_^
Name Size Date
060-1.fig270.39 kB01-08-12 10:26
060.bmp406.05 kB18-07-12 09:51
060.fig261.75 kB01-08-12 09:39
080.bmp406.05 kB18-07-12 09:51
090l.bmp406.05 kB18-07-12 09:51
090r.bmp406.05 kB18-07-12 09:51
090r.fig277.17 kB01-08-12 09:34
90.bmp406.05 kB20-07-12 16:14
circ.m533.00 B30-07-12 15:32
correct_image_061g.bmp3.82 MB01-08-12 19:52
fitellipse.m3.13 kB28-07-12 22:38
markpoint.m621.00 B01-08-12 09:57
rec.asv3.29 kB04-08-12 10:40
rec.m3.32 kB04-08-12 11:16
rec2.asv2.02 kB01-08-12 15:54
rec2.m2.02 kB01-08-12 15:55
shiyan5.m3.69 kB27-07-12 15:30
test.asv1.59 kB01-08-12 11:04
test.m1.63 kB01-08-12 16:43
test2.asv3.92 kB01-08-12 10:09
Thumbs.db36.00 kB18-07-12 09:51
trapezial_correction-wrong.m1.50 kB02-08-12 16:59
trapezial_correction.asv3.23 kB02-08-12 17:23
trapezial_correction.m1.66 kB04-08-12 10:09
Untitled.m1.45 kB26-07-12 20:43
Untitled3.m3.98 kB25-07-12 15:44
Untitled5.m3.75 kB23-07-12 08:58
Untitled6.m3.82 kB23-07-12 15:35
截图.jpg125.80 kB29-07-12 20:25
<校正>0.00 B09-04-13 09:07
...
Sponsored links

fitellipse.m (3.14 MB)

Need 2 point
Your Point(s)

Your Point isn't enough.

Get point immediately by PayPal

More(Debit card / Credit card / PayPal Credit / Online Banking)

Submit your source codes. Get more point

LOGIN

Don't have an account? Register now
Need any help?
Mail to: support@codeforge.com

切换到中文版?

CodeForge Chinese Version
CodeForge English Version

Where are you going?

^_^"Oops ...

Sorry!This guy is mysterious, its blog hasn't been opened, try another, please!
OK

Warm tip!

CodeForge to FavoriteFavorite by Ctrl+D