7
$\begingroup$

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.

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$
    – user484
    Commented Jul 21, 2014 at 9:06

5 Answers 5

Reset to default
10
$\begingroup$ 
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

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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 $\endgroup$
    – khanhnguyendata
    Commented Jun 27, 2014 at 5:50
  • $\begingroup$ @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. $\endgroup$
    – ubpdqn
    Commented Jun 27, 2014 at 6:00
  • $\begingroup$ @ubpdqn No apologies needed! If it weren't for you I wouldn't even know where to start to solve my problem. $\endgroup$
    – khanhnguyendata
    Commented Jun 27, 2014 at 6:31
  • $\begingroup$ * Larson 2011 Precalculus with Limits p. 761 $\endgroup$
    – khanhnguyendata
    Commented Jun 27, 2014 at 6:32
  • $\begingroup$ 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? $\endgroup$
    – khanhnguyendata
    Commented Jun 27, 2014 at 6:38
5
$\begingroup$

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];
Column[TraditionalForm /@ {ellipse, e2, e3}]

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

enter image description here

enter image description here

Improve this answer
$\endgroup$
5
$\begingroup$

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]

enter image description here

Improve this answer
$\endgroup$
5
$\begingroup$

Something very different (not sure when this function showed up):

ellipse = -1.43328 - 0.281261 x + 0.00228797 x^2 - 0.832702 y + 
   0.00138781 x y + 0.00209387 y^2;
ResourceFunctionHelpers`EllipseProperties[ellipse == 0, {x, y}]
(*
<|
 "Equation" -> -1.43328 - 0.281261 x + 0.00228797 x^2 - 0.832702 y + 
    0.00138781 x y + 0.00209387 y^2 == 0, 
 "Foci" -> {{-107.329, 323.282}, {109.906, 73.5497}}, 
 "Vertices" -> {{-154.738, 377.783}, {157.316, 19.0483}}, 
 "Covertices" -> {{-127.479, 86.4042}, {130.057, 310.427}}, 
 "Eccentricity" -> 0.696147,
 "Center" -> {1.28886, 198.416}, 
 "SemimajorAxisLength" -> 237.733,
 "SemiminorAxisLength" -> 170.668, 
 "AreaEnclosed" -> 127466.,
 "SemiLatusRectum" -> 122.523, 
 "Directrices" -> {
   y == 625.907 + 0.869872 x, 
   y == -56.5617 + 0.869872 x}, 
 "AxesOfSymmetry" -> {
   y == 199.897 - 1.14959 x, 
   y == 197.295 + 0.869872 x}, "Circumference" -> 1291.7, 
 "RotationAngle" -> -0.854878, 
 "BaseConic" -> <|
   "Equation" -> 0.0000176938 x^2 + 0.0000343316 y^2 == 1, 
   "RotationAngle" -> -0.854878, 
   "DisplacementFromOrigin" -> {1.28886, 198.416}|>|>
*)
Improve this answer
$\endgroup$
1
$\begingroup$

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

Mathematica graphics

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.