# An easier way to find the properties of an ellipse from its equation?

In this Q&A about fitting an ellipse to a set of points, there are multiple answers that generated general equations of the ellipse, like this one by @ubpdqn:

However, the steps to find out the properties of the ellipse (namely, major axis, minor axis, center, and rotation) from its given equation seems pretty complicated. Therefore, my goal is to try finding these properties in a more intuitive way (if possible).

Let's say I start with this general equation of a particular ellipse:

ellipse = -1.43328 - 0.281261 x + 0.00228797 x^2 - 0.832702 y +
0.00138781 x y + 0.00209387 y^2;
ContourPlot[ellipse == 0, {x, -500, 500}, {y, -500, 500}]


A way to find its properties is to start with a simple centered, axis-aligned ellipse with the general formula, where a and b represents the half-lengths of the ellipse' axes: $$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$$

The second step is to subtract the expression by 1 to convert our equation to the form expr=1 to expr=0, which is the way our general formula is written (see top of question).

Then, this ellipse is rotated a certain degree, which means that each x and y in the equation is corrected with the rotation angle alpha [Ref], like so: $$x\to x \cos (\alpha )-y \sin (\alpha ),y\to x \sin (\alpha )+y \cos (\alpha )$$

This tilted ellipse is then brought off-center to make our original ellipse, which will transform x and y like so: $$x\to x-h,y\to y-k$$ where (h,k) represents the coordinate of the ellipse' center.

This is the final formula as done via MMA:

ellipseTransform =

 Collect[#, {x,
y}] &@(x^2/a^2 + y^2/b^2 -
1 /. {x -> x Cos[\[Alpha]] - y Sin[\[Alpha]],
y -> x Sin[\[Alpha]] + y Cos[\[Alpha]]} /. {x -> x - h,
y -> y - k})


, which gives

A plot of this formula confirms that this is indeed an ellipse

ContourPlot[(ellipseTransform /. {a -> 300, b -> 200, h -> 100,
k -> 50, \[Alpha] -> 30 \[Degree]}) == 0, {x, -500,
500}, {y, -500, 500}, Axes -> True]


As you can see the transformed equation of the simple ellipse now resembles the form of our given ellipse, which includes a constant term, an x term, y term, xy term, x^2 term, and y^2 term (6 terms) in total. My goal now is to somehow solve for a, b, h, k and alpha using the fact that these 2 expressions are basically equivalent.

Now this is when I got stuck, as I've tried a few ways to equate these 2 expressions to look for what I want but to no avail:

Another (perhaps related) problem is that there are 5 unknowns (a,b,h,k,apha) but 6 equations i.e. the 6 coefficients in my ellipse equation. I suppose this problem is more a mathematics one than of Mathematica but I'm wondering if there's any way to solve (or at least optimize) an overdefined set of equations.

• I feel like the central question, "What is an intuitive way to find the major axis, minor axis, center, and rotation of an ellipse from its given quadratic equation?", is one better suited for the mathematics StackExchange. Once you have the mathematics, it should be short work to implement it in Mathematica. – Rahul Jul 21 '14 at 9:06

a = 0.00228797 ; b = 0.00138781/2; c = 0.00209387; d = -0.281261; e = -0.832702;
f = -1.43328;
ellipse1 = a x^2 + 2 b x y + c y^2 + d x + e y + f;
ContourPlot[ellipse1 == 0, {x, -500, 500}, {y, -500, 500}, ImageSize -> 200]


uv =First@Solve[ForAll[{x, y},
a (x - u)^2 + 2 b (x - u) (y - v) + c (y - v)^2 + f2 == ellipse1], {u, v, f2}, Reals]
ellipse2 = ellipse1 /. {x -> x + u, y -> y + v} /. uv // Simplify // Chop
ContourPlot[{ellipse1 == 0, ellipse2 == 0}, {x, -500, 500}, {y, -500, 500},
ImageSize -> 200, ContourStyle -> {Red, Blue}, Axes -> True, AxesOrigin -> {0, 0}]


solphi = First@Solve[Tan[2 phi] == 2 b/(a - c), phi] // Quiet
rotate = RotationMatrix[phi /. solphi].{x, y};
ellipse3 = ellipse2 /. {x -> First@rotate, y -> Last@rotate} // Simplify // Chop
ContourPlot[{ellipse1 == 0, ellipse2 == 0, ellipse3 == 0}, {x, -500, 500}, {y, -500,
500}, ImageSize -> 200, ContourStyle -> {Red, Blue, Green}, Axes -> True,
AxesOrigin -> {0, 0}]


Edit 1 Way 2

Clear["Global*"]
ellipse1 = -1.43328 - 0.281261 x + 0.00228797 x^2 - 0.832702 y +
0.00138781 x y + 0.00209387 y^2;
ellipse2 = Collect[#, {x, y}] &@(c (x^2/a^2 + y^2/b^2 - 1) /. {x ->
x Cos[\[Alpha]] - y Sin[\[Alpha]], y -> x Sin[\[Alpha]] + y Cos[\[Alpha]]} /.
{x -> x - h, y -> y - k});


It must has a “c” in this function,because if f1==0 and f2==0 then f1 == c f2 not f1==f2.

func = Flatten@MapThread[Equal, CoefficientList[#, {x, y}] & /@
{ellipse1, ellipse2}, 2]// DeleteCases[#, True] & // Simplify


There are six functions for six variables(a,b,c,h,k,[Alpha]).It is too difficult for NSolve to solve this six functions,so use FindRoot.

FindRoot[func, {a, 40}, {b, 50}, {c, 30}, {h, 2}, {k, 300}, {\[Alpha],Pi/4}]


• Brilliant! I was somewhat intimidated by your second step (rotation after centering) because matrices depress and confuse me, but I tried a more intuitive method using the transformation formula for rotation in the [link] (maa.org/external_archive/joma/Volume8/Kalman/General.html), also here* in my question and I got the same result: screengrab – seismatica Jun 27 '14 at 5:50
• @Cheminqi this is clearer. You use Mathematica's solve capabilities to the do the algebraic grunt work and then rotate the centered ellipse (I used eigenvectors). Apologies to seismatica for the opacity of my approach. – ubpdqn Jun 27 '14 at 6:00
• @ubpdqn No apologies needed! If it weren't for you I wouldn't even know where to start to solve my problem. – seismatica Jun 27 '14 at 6:31
• * Larson 2011 Precalculus with Limits p. 761 – seismatica Jun 27 '14 at 6:32
• Ok so I was curious if I could use the ForAll method from this answer directly to my original method: but Solve[ForAll[{x, y}, ellipseTransform == ellipse], {a, b, h, k, \[Alpha]}, Reals] has been running for an hour and it still hasn't returned an answer. Any tip on making it work? – seismatica Jun 27 '14 at 6:38

I post this for completeness (accepting that it was somewhat opaque).

In the series of questions it was known that the quadric form was an ellipse (this could be established using determinants as a previous poster has used).

Accepting an ellipse quadratic form,e.g. the test one provided

ellipse = -1.43328 - 0.281261 x + 0.00228797 x^2 - 0.832702 y +
0.00138781 x y + 0.00209387 y^2


The following:

cf[pol_] := Module[{c, h, k, cons},
c = CoefficientRules[pol];
{h, k} =
0.5 Inverse[{{-({2, 0} /. c), -0.5 ({1, 1} /.
c)}, {-0.5 ({1, 1} /. c), -({0, 2} /. c)}}].{{1, 0} /.
cr, {0, 1} /. c};
cons = {({0, 0} /.
c), ({2, 0} /. c) h^2, ({0, 2} /. c) k^2, ({1, 1} /.
c) h k}.{1, -1, -1, -1};
{h, k, cons};
Expand[({2, 0} /. c) (x)^2 + ({0, 2} /. c) (y)^2 + ({1, 1} /.
c) x y + cons]]
prax[pol_] := Module[{c, mat},
c = CoefficientRules[pol];
mat = {{{2, 0} /. c,
0.5 ({1, 1} /. c)}, {0.5 ({1, 1} /. c), {0, 2} /. c}};
Eigenvalues[mat].{x^2, y^2} + ({0, 0} /. c)
]


cf centers the ellipse. The ugly code is just completing the square and subtracting the constant to correct.

prax: just uses eigenvalues of matrix of quadric form to eliminate the cross term (x y).

Applying to test case:

e2 = cf[ellipse];
e3 = prax[e2];


This is consistent (reassuringly) with Chenminqis answer.

Visualizing:

ContourPlot[{ellipse == 0, e2 == 0, e3 == 0}, {x, -500,
500}, {y, -500, 500}, ContourStyle -> {Red, Green, Blue}]


Something little bit different.

a = 0.00228797; b =
0.00138781/
2; c = 0.00209387; d = -0.281261; e = -0.832702; f = -1.43328;
ellipse1 = a x^2 + 2 b x y + c y^2 + d x + e y + f;

point1 = {x, y} /. Maximize[{y, ellipse1 == 0}, {x, y}][[2]];
point2 = {x, y} /. Maximize[{-y, ellipse1 == 0}, {x, y}][[2]];
center = (point1 + point2)/2;
angle = phi /. Solve[Tan[2 phi] == 2 b/(a - c), phi][[1]] // Quiet;
ellipse2 = ellipse1 /. Thread[{x, y} -> {x, y} + center];
r = RotationTransform[angle];
ellipse3 = ellipse2 /. Thread[{x, y} -> r[{x, y}]];

ContourPlot[{ellipse1 == 0, ellipse2 == 0, ellipse3 == 0}, {x, -300,
300}, {y, -300, 500}, ImageSize -> 300, Axes -> True,
AspectRatio -> Automatic]


Noticing that for the general equation of an ellipse

$$a x^2+b y^2+c x y+d x+e y+f=0$$

we have $$c=0$$ if it's axis-aligned.

For the general equation of an axis-aligned ellipse $$a x^2+b y^2+d x+e y+f=0$$ we have $$d=e=0$$ if it's centered.

So this problem can be solved like this:

exp1= -1.43328 - 0.281261 x + 0.00228797 x^2 - 0.832702 y + 0.00138781 x y + 0.00209387 y^2;
newxy = RotationTransform[θ]@TranslationTransform[{a, b}]@{x, y};
exp2 = exp1 /. Thread[{x, y} -> newxy] // Expand;
angle = Solve[{Total@Cases[exp2, a_ x y -> a, Infinity] == 0, 0 < θ < Pi/2}, θ]
exp3 = exp2 /. angle // Chop // Expand;
vector = Solve[Total@Cases[exp3, #, Infinity] == 0 & /@ {a_ x -> a, b_ y -> b}]

ContourPlot[Evaluate[{0 == exp1, 0 == exp3 /. {a -> 0, b -> 0}, 0 == exp3 /. vector}],
{x, -500, 500}, {y, -500, 500}, Axes -> True]

{{θ -> 0.715919}}
{{a -> 131.195, b -> 148.857}}
`