Unnormalization of ellipse coefficients after direct ellipse fitting

Question

I am trying to understand the normalization and "unnormalization" steps in the direct least squares ellipse fitting algorithm developed by Fitzgibbon, Pilu and Fisher (improved by Halir and Flusser).

EDITED: More details about theory added. Is eigenvalue problem where confusion stems from?

Short theory:

The ellipse is represented by an implicit second order polynomial (general conic equation):

where:

To constrain this general conic to an ellipse, the coefficients must satisfy the quadratic constraint:

which is equivalent to:

where C is a matrix of zeros except:

The design matrix D is composed of all data points x sub i.

The minimization of the distance between a conic and the data points can be expressed by a generalized eigenvalue problem (some theory has been omitted):

Denoting:

We now have the system:

If we solve this system, the eigenvector corresponding to the single positive eigenvalue is the correct answer.

The code:

The code snippets here are directly from the MATLAB code provided by the authors: http://research.microsoft.com/en-us/um/people/awf/ellipse/fitellipse.html

The data input is a series of (x,y) points. The points are normalized by subtracting the mean and dividing by the standard deviation (in this case, computed as half the range). I'm assuming this normalization allows for a better fit of the data.

% normalize data
% X and Y are the vectors of data points, not normalized
mx = mean(X);
my = mean(Y);
sx = (max(X)-min(X))/2;
sy = (max(Y)-min(Y))/2; 

x = (X-mx)/sx;
y = (Y-my)/sy;

% Build design matrix
D = [ x.*x  x.*y  y.*y  x  y  ones(size(x)) ];

% Build scatter matrix
S = D'*D;   %'

% Build 6x6 constraint matrix
C(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;

[gevec, geval] = eig(S,C);

% Find the negative eigenvalue
I = find(real(diag(geval)) < 1e-8 & ~isinf(diag(geval)));

% Extract eigenvector corresponding to negative eigenvalue
A = real(gevec(:,I));

After this, the normalization is reversed on the coefficients:

par = [
  A(1)*sy*sy,   ...
  A(2)*sx*sy,   ...
  A(3)*sx*sx,   ...
  -2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy,   ...
  -A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy,   ...
  A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my   ...
  - A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my   ...
  + A(6)*sx*sx*sy*sy   ...
  ]';

At this point, I'm not sure what happened. Why is the unnormalization of the last three coefficients of A (d, e, f) dependent on the first three coefficients? How do you mathematically show where these unnormalization equations come from?

The 2 and 1 coefficients in the unnormalization lead me to believe the constraint matrix must be involved somehow.

Please let me know if more detail is needed on the method...it seems I'm missing how the normalization has propagated through the matrices and eigenvalue problem.

Any help is appreciated. Thanks!

Answer 1

At first, let me formalize the problem in a homogeneous space (as used in Richard Hartley and Andrew Zisserman's book Multiple View Geometry):
Assume that,

P=[X,Y,1]'

is our point in the unnormalized space, and

p=lambda*[x,y,1]'

is our point in the normalized space, where lambda is an unimportant free scale (in homogeneous space [x,y,1]=[10*x,10*y,10] and so on).

Now it is clear that we can write

x = (X-mx)/sx;
y = (Y-my)/sy;

as a simple matrix equation like:

p=H*P;  %(equation (1))

where

H=[1/sx,   0,     -mx/sx;
   0,      1/sy,  -my/sy;
   0,      0,          1];

Also we know that an ellipse with the equation

A(1)*x^2 + A(2)*xy + A(3)*y^2 + A(4)*x + A(5)*y + A(6) = 0  %(first representation)

can be written in matrix form as:

p'*C*p=0    %you can easily verify this by matrix multiplication

where

C=[A(1),    A(2)/2,  A(4)/2;
   A(2)/2,  A(3),    A(5)/2;
   A(4)/2,  A(5)/2,  A(6)];  %second representation

and

p=[x,y,1]

and it is clear that these two representations of an ellipse are exactly the same and equivalent.

Also we know that the vector A=[A(1),A(2),A(3),A(4),A(5),A(6)] is a type-1 representation of the ellipse in the normalized space.
So we can write:

p'*C*p=0

where p is the normalized point and C is as defined previously. Now we can use the "equation (1): p=HP" to derive some good result:

(H*P)'*C*(H*P)=0

=====>

P'*H'*C*H*P=0

=====>

P'*(H'*C*H)*P=0

=====>

P'*(C1)*P=0   %(equation (2))

We see that the equation (2) is an equation of an ellipse in the unnormalized space where C1 is the type-2 representation of ellipse and we know that:

C1=H'*C*H

Ans also, because the equation (2) is a zero equation we can multiply it by any non-zero number. So we multiply it by sx^2*sy^2 and we can write:

C1=sx^2*sy^2*H'*C*H

And finally we get the result

C1=[                                                          A(1)*sy^2,                                                         (A(2)*sx*sy)/2,                                                                                                                                                              (A(4)*sx*sy^2)/2 - A(1)*mx*sy^2 - (A(2)*my*sx*sy)/2;
                                                        (A(2)*sx*sy)/2,                                                              A(3)*sx^2,                                                                                                                                                              (A(5)*sx^2*sy)/2 - A(3)*my*sx^2 - (A(2)*mx*sx*sy)/2;
    -(- (A(4)*sx^2*sy^2)/2 + (A(2)*my*sx^2*sy)/2 + A(1)*mx*sx*sy^2)/sx,     -(- (A(5)*sx^2*sy^2)/2 + A(3)*my*sx^2*sy + (A(2)*mx*sx*sy^2)/2)/sy,     (mx*(- (A(4)*sx^2*sy^2)/2 + (A(2)*my*sx^2*sy)/2 + A(1)*mx*sx*sy^2))/sx + (my*(- (A(5)*sx^2*sy^2)/2 + A(3)*my*sx^2*sy + (A(2)*mx*sx*sy^2)/2))/sy + A(6)*sx^2*sy^2 - (A(4)*mx*sx*sy^2)/2 - (A(5)*my*sx^2*sy)/2]

which can be transformed into the type-2 ellipse and get the exact result we were looking for:

[ A(1)*sy^2, A(2)*sx*sy, A(3)*sx^2, A(4)*sx*sy^2 - 2*A(1)*mx*sy^2 - A(2)*my*sx*sy, A(5)*sx^2*sy - 2*A(3)*my*sx^2 - A(2)*mx*sx*sy,    A(2)*mx*my*sx*sy + A(1)*mx*my*sy^2 + A(3)*my^2*sx^2 + A(6)*sx^2*sy^2 - A(4)*mx*sx*sy^2 - A(5)*my*sx^2*sy]

If you are curious how I managed to caculate these time-consuming equations I can give you the matlab code to do it for you as follows:

syms sx sy mx my
syms a b c d e f
C=[a,   b/2, d/2;
   b/2, c,   e/2;
   d/2, e/2, f];    
H=[1/sx,  0,    -mx/sx;
   0,     1/sy, -my/sy;
   0,     0,     1];
C1=sx^2*sy^2*H.'*C*H
a=[Cp(1,1), 2*Cp(1,2), Cp(2,2), 2*Cp(1,3), 2*Cp(2,3), Cp(3,3)]