MATLAB: Ellipse implicit equation coefficients

algebracurve fittingellipsegeometryleast squares

Hi I have a set of 2D points and I want to know the coefficients of the following implicit equation of the ellipse: Ax^2+2Bxy+Cy^2=1

Any idea?

Moreover, I have found the best fitting ellipse with its major, minor axes, center and orientation. Can I can use this information to find the above mentioned three coefficients (A,B and C)?

The implementations of ellipse fitting already available in Matlab Central use the general form of ellipse, I do not need coefficients for those.

Best Regards Wajahat

Best Answer

First, if the centre is nonzero, then that form of the ellipse equation will not work. You need the more general quadratic form

(*) x^' * A * x + c' * x + d = 0

where A is a 2x2 symmetric matrix, c a 2x1 vector and d a scalar. Multiplying out gives the full equation

a_11 x^2 + 2a_12 xy + a_22 y^2 + c_1 x + c_2 y + d = 0

If you used my code: http://www.mathworks.com/matlabcentral/fileexchange/15125-fitellipse-m/ then what you get when you fit your ellipse is the centre z, the orientation of the major axis as a counterclockwise rotation from the x-axis alpha , and a and b the semimajor and semiminor axes respectively. Any of the other fitting codes probably return similar parameters.

These correspond to the following parametric equation for the ellipse (parametrised by theta in [0, 2\pi])

x = z + Q * [a * cos(theta); b * sin(theta)]

where

Q = [cos(alpha), -sin(alpha); sin(alpha) cos(alpha)]

All you need to do to get back to the quadratic form is to eliminate theta. First some algebra

diag([1/a, 1/b]) * Q'*(x - z) = [cos(theta); sin(theta)] = u

Using the fact that u'*u = 1,

(x - z)' * Q * diag([1/a^2, 1/b^2]) * Q' * (x - z) = 1

(x - z)' * A * (x - z) = 1

where A is symmetric. Expanding this out gives the required form:

x'*A*x - 2*z'*A*x + z'*A*z - 1 = 0

So linking back to my first equation (*), the MATLAB code for generating the matrices and vectors required are

Q = [cos(alpha), -sin(alpha); sin(alpha) cos(alpha)];
A = Q * diag([a^-2, b^-2]) * Q';
c = 2*A*z;
d = z'*A*z - 1;

Related Solutions

MATLAB: Fitting an ellipse to the binary iris image

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.
%
%See also ellipsoidfit
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;

MATLAB: Affine transformation that takes a given, known ellipse and maps it to a circle with diameter equal to the major axis.

I'm writing this down as a new answer so I can get markup for the code snippets. Let's work with the first blob from regionprops:

stats = regionprops(L)

'Orientation' is alpha, but in degrees:

alpha = pi/180 * stats(1).Orientation;
Q = [cos(alpha), -sin(alpha); sin(alpha), cos(alpha)];

'Centroid' is x0, but as a row vector (I think)

x0 = stats(1).Centroid.';

'MajorAxisLength' and 'MinorAxisLength' are a and b respectively

a = stats(1).MajorAxisLength;
b = stats(1).MinorAxisLength;

We've now got everything, so lets assemble these into the matrices we need:

S = diag([1, a/b]);
C = Q*S*Q';
d = (eye(2) - C)*x0;

You now have the matrices you need to build the affine transformation. Now if you're using maketform, it wants an transformation matrix that acts on homogeneous coordinates in row vector form, not columns, so we take the transpose of what you might expect.

tform = maketform('affine', [C d; 0 0 1]');

Then with a bit of luck, this will work (I'm doing this without the benefit of the image processing toolbox on my computer)

Im2 = imtransform(Im, tform)

Best Answer

Related Solutions

MATLAB: Fitting an ellipse to the binary iris image

MATLAB: Affine transformation that takes a given, known ellipse and maps it to a circle with diameter equal to the major axis.

Related Question