Standard form of ellipse

Question

I'm getting ellipses as level curves of a fit dataset. After selecting a particular ellipse, I would like to report it as a center point, semi-major and minor axes lengths, and a rotation angle. In other words, I would like to transform (using mathematica) my ellipse equation from the form:

Ax^2 + By^2 + Cx + Dy + Exy + F = 0

to a more standard form:

((xCos[alpha] - ySin[alpha] - h)^2)/(r^2) + ((xSin[alpha] + yCos[alpha] - k)^2)/(s^2) = 1

where (h,k) is the center, alpha is the rotation angle, and r and s are the semi-axes

The actual equation I'm attempting to transform is

 1.68052 x - 9.83173 x^2 + 4.89519 y - 1.19133 x y - 9.70891 y^2 + 6.09234 = 0

I know the center point is the fitted maximum, which is:

 {0.0704526, 0.247775}

Answer 1

I posted a version of this answer on Math SE since it benefits a lot from proper mathematical typesetting. The example there is simpler as well, and there are some extra details.

The following description follows the German Wikipedia article Hauptachsentransformation. Its English counterpart, according to inter-wiki links, is principal component analysis. I find the former article a lot more geometric than the latter. The latter has a strong focus on statistical data, though, so it might be useful for you nevertheless.

Rotation

Your ellipse is described as

        [A   E/2]   [x]             [x]
[x y] * [E/2   B] * [y]  +  [C D] * [y] + F = 0

First you identify the rotation. You do this by identifying the eigenvalues and eigenvectors of this 2×2 matrix. These eigenvectors will form an orthogonal matrix describing your rotation: its entries are the Sin[alpha] and Cos[alpha] from your formula.

With your numbers, you get

[A  E/2]   [-0.74248  0.66987]   [-10.369  0     ]   [-0.74248 -0.66987]
[E/2  B] = [-0.66987 -0.74248] * [ 0      -9.1715] * [ 0.66987 -0.74248]

The first of the three factors is the matrix formed by the eigenvectors, each normalized to unit length. The central matrix has the eigenvalues on the diagonal, and the last one is the transpose of the first. If you multiply the vector (x,y) with that last matrix, then you will change the coordinate system in such a way that the mixed term vanishes, i.e. the x and y axes are parallel to the main axes of your ellipse. This is just what happens in your desired formula, so now you know that

Cos[alpha] = -0.74248  (-0.742479398678 with more accuracy)
Sin[alpha] =  0.66987  ( 0.669868899516)

Translation

If you multiply the row vector [C D] in the above formula with the first of the three matrices, then this effect will exactly cancel the multiplication of (x, y) by the third matrix. Therefore in that changed coordinate system, you use the central diagonal matrix for the quadratic term, and this product for the linear term.

                     [-0.74248  0.66987]
[1.68052, 4.89519] * [-0.66987 -0.74248] = [-4.5269 -2.5089]

Now you have to complete the square independently for x and y, and you end up with a form from which you can read the center coordinates.

-10.369x² -4.5269x = -10.369(x + 0.21829)² + 0.49408
-9.1715y² -2.5089y = -9.1715(y + 0.13677)² + 0.17157

h = -0.21829  (-0.218286476695)
k = -0.13677  (-0.136774259156)

Note that h and k describe the center in the already rotated coordinate system; to obtain the original center you'd multiply again with the first matrix:

[-0.74248  0.66987]   [-0.21829]   [0.07045]
[-0.66987 -0.74248] * [-0.13677] = [0.24778]

which fits your description.

Scaling

The completed squares above contributed some more terms to the constant factor F:

6.09234 + 0.49408 + 0.17157 = 6.7580

Now you move this to the right side of the equation, then divide the whole equation by this number so that you get the = 1 from your desired form. Then you can deduce the radii.

1    -10.369
-- = ------- = 1.5344
r²   -6.7580

1    -9.1715
-- = ------- = 1.3571
s²   -6.7580

r = 0.80730  (0.807304599162099)
s = 0.85840  (0.858398019487315)

Verifying the result

Now let's check that we didn't make any mistakes. With the parameters we found, you can piece together the equation

  ((-0.74248*x - 0.66987*y + 0.21829)^2)/(0.80730^2)
+ (( 0.66987*x - 0.74248*y + 0.13677)^2)/(0.85840^2) = 1

Move the 1 to the left side, and multiply by -6.7580, and you should end up with the original equation. Expanding that (with the extra precision versions printed in parentheses), you'll get

-9.8317300000 x^2
-1.1913300000 x y
+1.6805200000 x
-9.7089100000 y^2
+4.8951900000 y
+6.0923400000

which is a perfect match for your input.

Answer 2

If you have h and k, you can use Lagrange Multipliers to maximize / minimize the function (x-h)^2+(y-k)^2 subject to the constraint of being on the ellipse. The maximum distance will be the major radius, the minimum distance the minor radius, and alpha will be how much they are rotated from horizontal.