Prev: Rounding error amplification in matrix/matrix multiplication
Next: MATLAB Ode45 needing help
From: Laura Suad on 9 Apr 2010 13:24
Hello,
I am using a ellipse fitting program. The one made by Tal Hendel, I need to have the: semimajor_axis, semiminor_axis, x0, y0, phi with dispersion (the fitting errors). How can I do to obtain these parameters?

From: Matt J on 9 Apr 2010 14:27
"Laura Suad" <lausuad(a)gmail.com> wrote in message <hpnnrk$34o$1(a)fred.mathworks.com>...
> Hello,
> I am using a ellipse fitting program. The one made by Tal Hendel, I need to have the: semimajor_axis, semiminor_axis, x0, y0, phi with dispersion (the fitting errors). How can I do to obtain these parameters?
===============

His/Her method doesn't seem to be a very good one. It ignores the noise in the matrix M. Below is a better way.

Aside from that, it is not clear what you are confused about. The Hendel code (and mine) give all the parameters that you requested as outputs, except for dispersion. For that, you'll need to specify how you would want to the fitting errors measured. One possibillity would be to find the perpendicular distances of your X-Y data to the ellipse...





function report=ellipsefit(XY)
%ELLIPSEFIT - form 2D ellipse fit to given x,y data
%
% report=ellipsefit(XY)
%
%in:
%
% XY: Input matrix of 2D coordinates to be fit. Each column XY(:,i) is [xi;yi]
%
%out: Finds the ellipse fitting the input data parametrized both as
% A*x^2+B*x*y C*y^2+D*x+E*y=1 and [x-x0,y-y0]*Q*[x-x0;y-y0]=1
%
% report: a structure output with the following fields
%
% report.Q: the matrix Q
% report.d: the vector [x0,y0]
% report.ABCDE: the vector [A,B,C,D,E]
% report.AxesDiams: The minor and major ellipse diameters
% report.theta: The counter-clockwise rotation of the ellipse.
%
%NOTE: The code will give errors if the data fit traces out a non-elliptic or
% degenerate conic section.


X=XY(1,:).';
Y=XY(2,:).';

M= [X.^2, X.*Y, Y.^2, X, Y, -ones(size(X,1),1)];

[U,S,V]=svd(M,0);
ABCDEF=V(:,end);

if size(ABCDEF,2)>1

error 'Data cannot be fit with unique ellipse'
else

ABCDEF=num2cell(ABCDEF);
end

[A,B,C,D,E,F]=deal(ABCDEF{:});


Q=[A, B/2;B/2 C];
x0=-Q\[D;E]/2;
dd=F+x0'*Q*x0;

Q=Q/dd;

[R,eigen]=eig(Q);
eigen=eigen([1,4]);

if ~all(eigen>=0), error 'Fit produced a non-elliptic conic section'; end

idx=eigen>0;
eigen(idx)=1./eigen(idx);
AxesDiams = 2*sqrt(eigen);

theta=atand(tand(-atan2(R(1),R(2))*180/pi));


report.Q=Q;
report.d=x0(:).';
report.ABCDE=[A, B, C, D, E]/F;
report.AxesDiams=sort(AxesDiams(:)).';
report.theta=theta;
 | 
Pages: 1
Prev: Rounding error amplification in matrix/matrix multiplication
Next: MATLAB Ode45 needing help