求助！如何用MAPLE求解对离散点的最优椭圆拟合

wb0330 · 发表于 2005-4-24 22:27:34

<

RE>如何用MAPLE求解对离散点的最优椭圆拟合？</PRE><

RE>请各位大侠告知！</PRE><

RE>这里有一个MATLAB的来自</PRE><PRE><a href="http://bbs.dartmouth.edu/~fangq/wiki/index.cgi?MathTools_FAQ" target="_blank" >http://bbs.dartmouth.edu/~fangq/wiki/index.cgi?MathTools_FAQ</A></PRE><PRE>function a = fitellipse(X,Y)</PRE><PRE> </PRE><PRE>% FITELLIPSE Least-squares fit of ellipse to 2D points.</PRE><PRE>% A = FITELLIPSE(X,Y) returns the parameters of the best-fit</PRE><PRE>% ellipse to 2D points (X,Y).</PRE><PRE>% The returned vector A contains the center, radii, and orientation</PRE><PRE>% of the ellipse, stored as (Cx, Cy, Rx, Ry, theta_radians)</PRE><PRE>%</PRE><PRE>% Authors: Andrew Fitzgibbon, Maurizio Pilu, Bob Fisher</PRE><PRE>% Reference: "Direct Least Squares Fitting of Ellipses", IEEE T-PAMI, 1999</PRE><PRE>%</PRE><PRE>% This is a more bulletproof version than that in the paper, incorporating</PRE><PRE>% scaling to reduce roundoff error, correction of behaviour when the input</PRE><PRE>% data are on a perfect hyperbola, and returns the geometric parameters</PRE><PRE>% of the ellipse, rather than the coefficients of the quadratic form.</PRE><PRE>%</PRE><PRE>% Example: Run fitellipse without any arguments to get a demo</PRE><PRE>if nargin == 0</PRE><PRE> % Create an ellipse</PRE><PRE> t = linspace(0,2);</PRE><PRE> </PRE><PRE> Rx = 300</PRE><PRE> Ry = 200</PRE><PRE> Cx = 250</PRE><PRE> Cy = 150</PRE><PRE> Rotation = .4 % Radians</PRE><PRE> </PRE><PRE> x = Rx * cos(t);</PRE><PRE> y = Ry * sin(t);</PRE><PRE> nx = x*cos(Rotation)-y*sin(Rotation) + Cx;</PRE><PRE> ny = x*sin(Rotation)+y*cos(Rotation) + Cy;</PRE><PRE> % Draw it</PRE><PRE> plot(nx,ny,'o');</PRE><PRE> % Fit it</PRE><PRE> fitellipse(nx,ny)</PRE><PRE> % Note it returns (Rotation - pi/2) and swapped radii, this is fine.</PRE><PRE> return</PRE><PRE>end</PRE><PRE> </PRE><PRE>% normalize data</PRE><PRE>mx = mean(X);</PRE><PRE>my = mean(Y);</PRE><PRE>sx = (max(X)-min(X))/2;</PRE><PRE>sy = (max(Y)-min(Y))/2;</PRE><PRE> </PRE><PRE>x = (X-mx)/sx;</PRE><PRE>y = (Y-my)/sy;</PRE><PRE> </PRE><PRE>% Force to column vectors</PRE><PRE>x = x(

;</PRE><PRE>y = y(

;</PRE><PRE> </PRE><PRE>% Build design matrix</PRE><PRE>D = [ x.*x x.*y y.*y x y ones(size(x)) ];</PRE><PRE> </PRE><PRE>% Build scatter matrix</PRE><PRE>S = D'*D;</PRE><PRE> </PRE><PRE>% Build 6x6 constraint matrix</PRE><PRE>C(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;</PRE><PRE> </PRE><PRE>% Solve eigensystem</PRE><PRE>[gevec, geval] = eig(S,C);</PRE><PRE> </PRE><PRE>% Find the negative eigenvalue</PRE><PRE>I = find(real(diag(geval)) < 1e-8 & ~isinf(diag(geval)));</PRE><PRE> </PRE><PRE>% Extract eigenvector corresponding to negative eigenvalue</PRE><PRE>A = real(gevec(:,I));</PRE><PRE> </PRE><PRE>% unnormalize</PRE><PRE>par = [</PRE><PRE> A(1)*sy*sy, ...</PRE><PRE> A(2)*sx*sy, ...</PRE><PRE> A(3)*sx*sx, ...</PRE><PRE> -2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy, ...</PRE><PRE> -A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy, ...</PRE><PRE> A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my ...</PRE><PRE> - A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my ...</PRE><PRE> + A(6)*sx*sx*sy*sy ...</PRE><PRE> ]';</PRE><PRE> </PRE><PRE>% Convert to geometric radii, and centers</PRE><PRE> </PRE><PRE>thetarad = 0.5*atan2(par(2),par(1) - par(3));</PRE><PRE>cost = cos(thetarad);</PRE><PRE>sint = sin(thetarad);</PRE><PRE>sin_squared = sint.*sint;</PRE><PRE>cos_squared = cost.*cost;</PRE><PRE>cos_sin = sint .* cost;</PRE><PRE> </PRE><PRE>Ao = par(6);</PRE><PRE>Au = par(4) .* cost + par(5) .* sint;</PRE><PRE>Av = - par(4) .* sint + par(5) .* cost;</PRE><PRE>Auu = par(1) .* cos_squared + par(3) .* sin_squared + par(2) .* cos_sin;</PRE><PRE>Avv = par(1) .* sin_squared + par(3) .* cos_squared - par(2) .* cos_sin;</PRE><PRE> </PRE><PRE>% ROTATED = [Ao Au Av Auu Avv]</PRE><PRE> </PRE><PRE>tuCentre = - Au./(2.*Auu);</PRE><PRE>tvCentre = - Av./(2.*Avv);</PRE><PRE>wCentre = Ao - Auu.*tuCentre.*tuCentre - Avv.*tvCentre.*tvCentre;</PRE><PRE> </PRE><PRE>uCentre = tuCentre .* cost - tvCentre .* sint;</PRE><PRE>vCentre = tuCentre .* sint + tvCentre .* cost;</PRE><PRE> </PRE><PRE>Ru = -wCentre./Auu;</PRE><PRE>Rv = -wCentre./Avv;</PRE><PRE> </PRE><PRE>Ru = sqrt(abs(Ru)).*sign(Ru);</PRE><PRE>Rv = sqrt(abs(Rv)).*sign(Rv);</PRE><PRE> </PRE><PRE>a = [uCentre, vCentre, Ru, Rv, thetarad];</PRE>

		自动登录	找回密码
密码			注-册-帐-号