7

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:

enter image description here

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}]

Mathematica graphics

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

Mathematica graphics

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]

Mathematica graphics

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:

Mathematica graphics

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.

Community
  • 1
seismatica
  • 5,101
  • 1
  • 22
  • 33
  • 1
    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. –  Jul 21 '14 at 09:06

5 Answers5

10
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]

enter image description here

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}]

enter image description here

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}]

enter image description here

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

enter image description here

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}]

enter image description here

Apple
  • 3,663
  • 16
  • 25
  • 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] (http://www.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 05: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 06: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 06:31
  • Larson 2011 Precalculus with Limits p. 761
    • – seismatica Jun 27 '14 at 06: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 06:38
  • Thank you again for the brilliant edit. Now I know exactly where I went wrong and how I could have done it. I applied your coefficient comparison method to my original equation--added a c term of course--and I got the same answer. – seismatica Jun 27 '14 at 16:34