Matlab routine for ellipse fitting

AaronSpark

发信人: Mars (FangQ), 信区: MathTools<BR>标  题: Matlab routine for ellipse fitting<BR>发信站: 达摩BigGreen BBS (Wed Sep 25 13:08:29 2002), 站内信件<BR><BR>From: tom (tom2959@21cn.com)<BR>Subject: ellipse fitting? <BR>Newsgroups: comp.lang.idl-pvwave<BR>View: Complete Thread (4 articles) | Original Format<BR>Date: 2002-04-27 06:42:17 PST <BR><BR><BR>I found a matlab function for ellipse, but it is not easy for me translate<BR>to IDl. For example,<BR><BR> % Solve eigensystem<BR>    [gevec, geval] = eig(S,C);<BR><BR>are there any function like eig(S,C) in IDL?<BR><BR>The matlab for ellips fitting is as following, who have a idl version?<BR><BR><BR><BR>function a = fitellipse(X,Y)<BR><BR>% FITELLIPSE  Least-squares fit of ellipse to 2D points.<BR>%        A = FITELLIPSE(X,Y) returns the parameters of the best-fit<BR>%        ellipse to 2D points (X,Y).<BR>%        The returned vector A contains the center, radii, and orientation<BR>%        of the ellipse, stored as (Cx, Cy, Rx, Ry, theta_radians)<BR>%<BR>% Authors: Andrew Fitzgibbon, Maurizio Pilu, Bob Fisher<BR>% Reference: "Direct Least Squares Fitting of Ellipses", IEEE T-PAMI, 1999<BR>%<BR>% This is a more bulletproof version than that in the paper, incorporating<BR>% scaling to reduce roundoff error, correction of behaviour when the input<BR>% data are on a perfect hyperbola, and returns the geometric parameters<BR>% of the ellipse, rather than the coefficients of the quadratic form.<BR>%<BR>%  Example:  Run fitellipse without any arguments to get a demo<BR>if nargin == 0<BR>  % Create an ellipse<BR>  t = linspace(0,2);<BR><BR>  Rx = 300<BR>  Ry = 200<BR>  Cx = 250<BR>  Cy = 150<BR>  Rotation = .4 % Radians<BR><BR>  x = Rx * cos(t);<BR>  y = Ry * sin(t);<BR>  nx = x*cos(Rotation)-y*sin(Rotation) + Cx;<BR>  ny = x*sin(Rotation)+y*cos(Rotation) + Cy;<BR>  % Draw it<BR>  plot(nx,ny,'o');<BR>  % Fit it<BR>  fitellipse(nx,ny)<BR>  % Note it returns (Rotation - pi/2) and swapped radii, this is fine.<BR>  return<BR>end<BR><BR>% normalize data<BR>mx = mean(X);<BR>my = mean(Y);<BR>sx = (max(X)-min(X))/2;<BR>sy = (max(Y)-min(Y))/2;<BR><BR>x = (X-mx)/sx;<BR>y = (Y-my)/sy;<BR><BR>% Force to column vectors<BR>x = x(:);<BR>y = y(:);<BR><BR>% Build design matrix<BR>D = [ x.*x  x.*y  y.*y  x  y  ones(size(x)) ];<BR><BR>% Build scatter matrix<BR>S = D'*D;<BR><BR>% Build 6x6 constraint matrix<BR>C(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;<BR><BR>% Solve eigensystem<BR>[gevec, geval] = eig(S,C);<BR><BR>% Find the negative eigenvalue<BR>I = find(real(diag(geval)) &lt; 1e-8 &amp; ~isinf(diag(geval)));<BR><BR>% Extract eigenvector corresponding to negative eigenvalue<BR>A = real(gevec(:,I));<BR><BR>% unnormalize<BR>par = [<BR>  A(1)*sy*sy,   ...<BR>      A(2)*sx*sy,   ...<BR>      A(3)*sx*sx,   ...<BR>      -2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy,   ...<BR>      -A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy,   ...<BR>      A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my   ...<BR>      - A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my   ...<BR>      + A(6)*sx*sx*sy*sy   ...<BR>      ]';<BR><BR>% Convert to geometric radii, and centers<BR><BR>thetarad = 0.5*atan2(par(2),par(1) - par(3));<BR>cost = cos(thetarad);<BR>sint = sin(thetarad);<BR>sin_squared = sint.*sint;<BR>cos_squared = cost.*cost;<BR>cos_sin = sint .* cost;<BR><BR>Ao = par(6);<BR>Au =   par(4) .* cost + par(5) .* sint;<BR>Av = - par(4) .* sint + par(5) .* cost;<BR>Auu = par(1) .* cos_squared + par(3) .* sin_squared + par(2) .* cos_sin;<BR>Avv = par(1) .* sin_squared + par(3) .* cos_squared - par(2) .* cos_sin;<BR><BR>% ROTATED = [Ao Au Av Auu Avv]<BR><BR>tuCentre = - Au./(2.*Auu);<BR>tvCentre = - Av./(2.*Avv);<BR>wCentre = Ao - Auu.*tuCentre.*tuCentre - Avv.*tvCentre.*tvCentre;<BR><BR>uCentre = tuCentre .* cost - tvCentre .* sint;<BR>vCentre = tuCentre .* sint + tvCentre .* cost;<BR><BR>Ru = -wCentre./Auu;<BR>Rv = -wCentre./Avv;<BR><BR>Ru = sqrt(abs(Ru)).*sign(Ru);<BR>Rv = sqrt(abs(Rv)).*sign(Rv);<BR><BR>a = [uCentre, vCentre, Ru, Rv, thetarad];<BR><BR><BR><BR>Google Home - Advertise with Us - Search Solutions - Language Tools - Jobs, Press, &amp; Help<BR><BR>?2002 Google<BR><BR>