Matlab routine for ellipse fitting
发信人: 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> = 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>% FITELLIPSELeast-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.*xx.*yy.*yxyones(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> = eig(S,C);<BR><BR>% Find the negative eigenvalue<BR>I = find(real(diag(geval)) < 1e-8 & ~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 = <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 = ;<BR><BR><BR><BR>Google Home - Advertise with Us - Search Solutions - Language Tools - Jobs, Press, & Help<BR><BR>?2002 Google<BR><BR>
页:
[1]