Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have 24 x-y points that are supposed to form an ellipse.

points={{31.4799,432.849},{-195.826,356.419},{210.029,121.779},{76.1586,-4.47992},{26.6061,-6.97711},{160.236,27.0398},{203.248,241.865},{-225.351,218.912},{-159.899,109.93},{-106.465,58.164},{-1.98952*10^-11,442.},{-229.146,324.489},{-31.4799,9.15114},{195.826,85.5807},{-210.029,320.221},{-76.1586,446.48},{-26.6061,448.977},{-160.236,414.96},{-203.248,200.135},{225.351,223.088},{159.899,332.07},{106.465,383.836},{-1.98952*10^-11,1.83076*10^-11},{229.146,117.511}};
points//ListPlot

points

I want to numerically fit an ellipse around those points. I tried to use the NMinimize method with least-square shown in this Q&A but the issue in this problem is that x and y are related parametrically (through t), and not directly like that other problem was. Worse, the ellipse will be off-center and tilted, so I can't use the nice $$\left(\frac{x}{a}\right)^2+\left(\frac{y}{b}\right)^2=1$$ but had to use this equation I found on Wikipedia (surprisingly without citation!):

ellipsewiki

My strategy was to generate 24 X(t) equations and 24 Y(t) equations, or more specifically, X(t1), Y(t1), ..., X(t24), X(t24) with xc, yc, a, b in them. Then, I would use NMinimize to optimize all of the square differences between my points and the ellipse. I would eventually get back xc, yc, a, b, CurlyPhi, and every single t values from t1 to t24. I tried to find a more elegant way of using NMinimized for parametric equations, and I believe the multivariate examples of NMinimized in the documentation are not exactly parametric equations. See also the linked Q&A above for an example.

Here's my code to find xc, yc, a, b, t1, ... , t24:

tList = ToExpression["t" <> ToString[#]] & /@ Range[1, 24];
equations = Join[
   MapThread[#1 - xc - a Cos[#2] Cos[\[CurlyPhi]] + 
      b Sin[#2] Sin[\[CurlyPhi]] &, {points[[All, 1]], tList}],
   MapThread[#1 - yc - a Cos[#2] Sin[\[CurlyPhi]] - 
      b Sin[#2] Cos[\[CurlyPhi]] &, {points[[All, 2]], tList}]];
solution = NMinimize[
  #.# &[equations],
  Join[tList, {xc, yc, a, b, \[CurlyPhi]}]]

equations

And this is what the ellipse looks like:

Show[

 (* Tilted ellipse/disk drawn using Disk and Rotate *)
 Graphics[
  {Yellow, 
   Rotate[Disk[{xc, yc}, {Abs@a, b}], \[CurlyPhi]] /. 
    solution[[2, 25 ;; 29]]}],

 (* Starting points *)
 points // ListPlot,

 (* Ellipse outline plotted by the 2 parametric equations *)
 ParametricPlot[{xc + a Cos[t] Cos[\[CurlyPhi]] - 
     b Sin[t] Sin[\[CurlyPhi]], 
    yc + a Cos[t] Sin[\[CurlyPhi]] + b Sin[t] Cos[\[CurlyPhi]]} /. 
   solution[[2, 25 ;; 29]], {t, 0, 2 \[Pi]}, PlotStyle -> Red],

 (* Point on ellipse corresponding to t1 to t24 *)
 ListPlot[{xc + a Cos[#] Cos[\[CurlyPhi]] - b Sin[#] Sin[\[CurlyPhi]],
       yc + a Cos[#] Sin[\[CurlyPhi]] + b Sin[#] Cos[\[CurlyPhi]]} /. 
     solution[[2, 25 ;; 29]] & /@ solution[[2, 1 ;; 24, 2]], 
  PlotStyle -> Red],

 Axes -> True]

ellipse

Even though I got what I want, I just feel that this is clunky way to solve the problem, not least because I don't need to solve for t1 to t24--I just needed xc, yc, a, and b. Another issue is that for some reason the value of my a (supposed to be the long axis) is negative and smaller in magnitude than the other axis*, and the rotate angle of the ellipse is crazy high (88.7 radian). I tried to apply some constrains to NMinimize (a>0 or 90 [Degree] < [CurlyPhi] < 150 [Degree]) but kept getting errors.

NMinimized error

*

ellipse wiki reference

Lastly, I'd really appreciate any general comment to make my code better as I'm still a beginner in MMA (you probably could tell by my zeal of pure functions in this example; I've just been reading a chapter about them).

share|improve this question
2  
This might be interesting for you: link –  eldo Jun 25 at 11:28
1  
i did this some years ago as a two level fit algorithm. Estimate xc,yc. Transform data to polar coordinates and use NonlinerarModelFit to fit the three parameters a,b,phi to the polar ellipse equation. Wrap that in a fit to optimize the center position. (my experience this worked far better than a direct 5 parameter fit) –  george2079 Jun 25 at 12:44
2  
(1) There are a number of approaches posted that seem good. About the only thing to be on the alert for is a "solution" that actually gives a hyperbola. This can happen in extreme cases of troublesome data. So you may want to enforce an inequality constraint in any optimization. –  Daniel Lichtblau Jun 25 at 18:41
3  
(2) Also maybe have a look at "Direct and specific least-square fitting of hyperbolae and ellipses" by P. O'Leary and P. Zsomber-Murray, Journal of Electronic Imaging 13(3):492-503, 2004. –  Daniel Lichtblau Jun 25 at 18:42
1  
See also the well-written description of ellipse parameters at mathworld.wolfram.com/Ellipse.html –  John McGee Jul 2 at 10:47

5 Answers 5

up vote 19 down vote accepted

Here is another approach. It could be improved (I am sure) to properly determined the principle axes and translation (if I get time I will aim to update):

lin = {#1^2, #1, #2, 2 #1 #2, #2^2} & @@@ points;
lm = LinearModelFit[lin, {1, a, b, c, d}, {a, b, c, d}]

Exploring model:

lm["ParameterTable"]

enter image description here

Determining quadric formula:

pa = lm["BestFitParameters"];
w[x_, y_] := pa.{1, x^2, x, y, 2 x y} - y^2;
ContourPlot[w[x, y] == 0, {x, -400, 400}, {y, -200, 500}, 
 Epilog -> Point[points]]

enter image description here

Formula:

TraditionalForm[-w[x, y] == 0]

enter image description here

UPDATE

I post this update to determine the translation from origin, the principle axes and the $a$ and $b$ of the desired form $\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$ and hence assessment of area of ellipse: viz $\pi a b$ etc. The axes can be swapped as desired.

The translation can be derived as follows:

{dx, dy} = 
 0.5 Inverse[{{-pa[[2]], -pa[[5]]}, {-pa[[5]], 1}}].{pa[[3]], pa[[4]]}

yielding:

{0.0000180983, 221.}

The axes can be derived from the matrix of quadric form:

mat = {{-pa[[2]], -pa[[5]]}, {-pa[[5]], 1}};
const = (-pa[[2]] dx^2 + dy^2 - 2 pa[[5]] dx dy) - pa[[1]]
trf[x_, y_] := -pa[[2]] (x - dx)^2 - 
  2 pa[[5]] (x - dx) (y - dy) + (y - dy)^2 - const

Then translating ellipse back to origin:

tran = Simplify@trf[x + dx, y + dy]

yields:

-49550.5 + 1.02504 x^2 + 0.498521 x y + y^2=0

Then looking at the eigensystem of matrix will allow visualization of principal axes:

fun[poly_, a_] := Module[{mat, val, vc, vcn, gr},
   mat = {{#1, #2/2}, {#2/2, #3}} & @@ (Coefficient[
       poly, {x^2, x y, y^2}]);
   {val, vc} = Eigensystem[mat];
   vcn = Normalize /@ vc;
   gr = Graphics[{Line[{{0, 0}, Sqrt[a] vcn[[1]]/Sqrt[Abs@val[[1]]]}],
       Line[{{0, 0}, Sqrt[a] vcn[[2]]/Sqrt[Abs@val[[2]]]}]}];
   Show[ContourPlot[poly == a, {x, -500, 500}, {y, -500, 500}], gr]
   ];

The values of $a$ and $b$ can e derived:

{ev, vec} = Eigensystem[mat];
st = {x^2, y^2}.ev;
{a, b} = Sqrt[1/Coefficient[Expand[st/const], {x^2, y^2}]]

yields:

{198.143, 254.846}

Putting visualizations all together:

cntr = ContourPlot[trf[x, y] == 0, {x, -500, 500}, {y, -200, 500}, 
  Epilog -> {Point[points], {Red, PointSize[0.03], Point[{dx, dy}]}}];
axes = fun[1.0250354570455185` x^2 + 0.49852055996040495` x y + y^2, 
  49550.46446156194`];
norm = ContourPlot[st == const, {x, -500, 500}, {y, -500, 500}];

These graphics and corresponding equations were threaded then exported as an animated gif:

grap = Show[#, PlotRange -> {{-500, 500}, {-500, 500}}] & /@ {cntr, 
   axes, norm}
eqns = Style[#, 20] & /@ {TraditionalForm[trf[x, y] == 0], 
   TraditionalForm[
    1.0250354570455185` x^2 + 0.49852055996040495` x y + y^2 == 
     49550.46446156194`], TraditionalForm[ell[x, y] == 0]}

enter image description here

The gif cycles through (i) data with fit and red point is {dx,dy} (ii) the second is the ellipse translated to the origin (iii) the ellipse axes are 'aligned' to cartesian axes.

share|improve this answer
    
Beautifully terse. You could make it more so with {#^2, ##, 2 ##, #2^2} & @@@ points. :-) –  Mr.Wizard Jun 25 at 16:09
    
@Mr.Wizard thank you for kind comment, edit and education –  ubpdqn Jun 26 at 9:30
    
Could you provide a reference for {dx, dy} = 0.5 Inverse[{{-pa[2], -pa[[5]]}, {-pa[[5]], 1}}].{pa[3], pa[4]}? I'm struggling to find a way to convert between general formula of an ellipse (x^2,y^2,xy,...) to the form that shows a and b. –  seismatica Jun 26 at 19:26
    
@seismatica {dx,dy} is just equating coefficients: $ax^2+bx+2cxy++dy+y^2= a(x-dx)^2+2c(x-dx)(y-dy)+(y-dy)^2-a dx^2-dy^2-2dx dy$ just solving using linear algebra. Daniel Lichtblau;s answer is a professional general way cf my hamfisted approach. –  ubpdqn Jun 26 at 22:50

In analytic geometry, the ellipse is defined as the set of points (X,Y) of the Cartesian plane that, in non-degenerate cases, satisfy the implicit equation enter image description here

with enter image description here and enter image description here where

enter image description here

enter image description here

enter image description here

Lets fit points with second-order curve (which include ellipse).

elipse = a11*x^2 + a22*y^2 + 2*a12*x*y + 2*a13*x + 2*a23*y + a33;
coeff = {a11, a22, a12, a13, a23, a33};
fitResult=FindFit[points/.{x_,y_}:>{x,y,0},{elipse,{#>0||#<0}&/@coeff},coeff,{x,y}];
elipse=elipse/.fitResult
(*{a11->0.002852802341,a22->0.002611181564,a12->0.0008537122851,a13->-0.188169993,a23->-0.5733778586,a33->-2.326357837}*)

Then check result

Show[{ContourPlot[res==0,{x,Min[points],Max[points]},{y,Min[points],Max[points]}],ListPlot[points]},PlotRange->All]

enter image description here

share|improve this answer
    
It becomes easier, and stabler, if you get rid of a23 and instead force the constant to be 1. –  Daniel Lichtblau Jun 26 at 0:06

Here is a standard direct way to get the principal exes and other transformation data. Find the mean of the points, subtract it to center them, and take the singular value decomposition. The third and second components thereof give the rotation and scaling data necessary to form a circle on which the first component, viewed as a point set, roughly lies. The average distance from origin of this new point set gives a radius. Now invert all that to get the ellipse equation.

While I'm sure that made no sense (least of all to me), the coding is straightforward.

mean = Mean[points];
newpts = Map[# - mean &, points];
{uu, ww, vv} = SingularValueDecomposition[newpts, 2]

(* Out[195]= {{{0.154787303139, 
   0.24681671674}, {0.264618646868, -0.0802782275335}, \
{-0.244940409458, 
   0.132178292936}, {-0.2488804106, -0.141810231081}, \
{-0.21297321559, -0.198934383898}, {-0.286761033433, \
-0.0175688333828}, {-0.138307266939, 
   0.244940128897}, {0.171124630754, -0.250504589808}, \
{0.0288583749054, -0.287472943299}, {-0.0558601605337, \
-0.280401389223}, {0.186667941, 
   0.221277578296}, {0.263211882778, -0.148978381541}, \
{-0.154787157436, -0.246816544024}, {-0.264618872812, 
   0.0802779596975}, {0.244940436909, -0.132178260395}, \
{0.248880505624, 0.141810343723}, {0.21297315013, 
   0.1989343063}, {0.286760891953, 
   0.0175686656724}, {0.13830729439, -0.244940096356}, \
{-0.171124603303, 0.250504622349}, {-0.0288583474542, 
   0.28747297584}, {0.0558601879848, 
   0.280401421764}, {-0.186667913548, -0.221277545756}, \
{-0.263211855327, 0.148978414082}}, {{876.375887309, 0.}, {0., 
   671.508170856}}, {{-0.672351542658, 
   0.740231992746}, {0.740231992746, 0.672351542658}}} *)

Here is the circle, centered at the origin, that the new points (from the left factor, uu), occupy.

ListPlot[uu, AspectRatio -> Automatic]

enter image description here

Now get the radius squared of this new circle.

rsqr = Mean[Map[#.# &, uu]]

(* Out[191]= 0.0833333333333 *)

To get to this origin-centered circle we run the transformations on our variables, {x,y}, to get new coordinates in terms thereof.

{nx, ny} = Inverse[vv.ww].({x, y} - mean);

So here is the ellipse equation.

expr = Expand[nx^2 + ny^2] == rsqr

(* Out[201]= 
0.0838086622042 - 0.000201425757605 x + 1.80374774713*10^-6 x^2 - 
  0.00075844948748 y + 9.11428834461*10^-7 x y + 
  1.71594919287*10^-6 y^2 == 0.0833333333333 *)

We'll check this pictorially.

ContourPlot[Evaluate@expr, {x, -400, 400}, {y, -200, 500}, 
 Epilog -> {Red, PointSize[Medium], Point[points]}]

enter image description here

share|improve this answer
    
I certainly enjoyed reading your code. Elegant and understandable. –  eldo Jun 25 at 20:26
    
@DanielLichtblau thank you for such a clear exposition of the linear algebra...closer and general cf contortions (which really use linear algebra in a tortured way) –  ubpdqn Jun 26 at 21:49

The approach mentioned in my comment using a polar representation:

 r[ a_, b_, theta_, t_] := 
     a b/Sqrt[b^2 Cos[theta] Sin[2 t] - a^2 Sin[2 t] Sin[theta]^2 +
          (a^2 + b^2) Sin[t]^2 Sin[theta]^2 +
          Cos[t]^2 (b^2 Cot[theta]^2 + a^2 Sin[theta]^2)];

some random sample data:

 data = With[ {a = 1, b = 1/2, theta = Pi/3 , x0 = 1, y0 = 1} ,
     Table[{x0, y0} + {Cos[t], Sin[t] } r[a, b, theta, t] +
     RandomVariate[NormalDistribution[0, .05], 2] ,
       {t, RandomVariate[UniformDistribution[{0, 2 Pi}], 100]} ]];

polar coordinate fit about a given center:

bestpfit[{xc_?NumericQ, yc_?NumericQ}, data_] := Module[{polardata,fit},
  polardata = { ArcTan @@ # , Norm[#] } & /@ ( # - {xc, yc} & /@ data);
  fit = NonlinearModelFit[ polardata, 
        r[ a, b, theta, t] , {a, b, theta } , { t} ];
  param = fit["BestFitParameters"];
  Norm@fit["FitResiduals"] ]

it helps to make a good guess at the center..then find xo,y0 to minimize the residual of the polar fit:

 guess = Mean[data];
 min = FindMinimum[ 
   bestpfit[{x0, y0}, data] , {{x0, guess[[1]]}, {y0, guess[[2]]}}] ;
 result = param~Join~min[[2]]

( note param is a global set by the last call to bestpfit .. probably there is a better way to do that )

{a -> 0.988222, b -> 0.475728, theta -> 0.990139, x0 -> 1.01157, y0 -> 1.00329}

(recall input was {1,1/2,1.0472,1,1})

 Show[{ListPlot[data, PlotStyle -> Red], 
      ParametricPlot[ {x0, y0} + {Cos[t], Sin[t] } r[a, b, theta, t] /. 
        result , {t, 0, 2 Pi}]}, PlotRange -> All]

enter image description here

This works well even with data not fully populating the ellipese:

 data = Select[ With[ {a = 1, b = 1/10, theta = Pi/6 , x0 = 1, y0 = 1} ,
      Table[{x0, y0} + {Cos[t], Sin[t] } r[a, b, theta, t] +
         RandomVariate[NormalDistribution[0, .01], 2] ,
           {t, RandomVariate[UniformDistribution[{0, 2 Pi}], 
                200]} ]] , #[[1]] > 1 && #[[2]] > 1/2 & ];

enter image description here

share|improve this answer
    
You've got my vote. I'm not sure why nobody else seems to have noticed this answer. –  Mr.Wizard Jun 25 at 17:57
    
You said in your comment that "this worked far better than a direct 5 parameter fit" in your experience. Can you elaborate a little on why that might be, or in what cases would this approach work better than the simpler direct fit? –  Rahul Narain Jun 25 at 18:45
    
It may well be that my 5 parameter algorithm wasn't very good (!). That was long ago and w/o mathematica. I'm sure the issue was fit success, not performance though. If I had time I'd try the other answers here on my last poor-data example.. –  george2079 Jun 25 at 19:10
    
+1 like this approach also –  ubpdqn Jun 25 at 22:23
Needs["MultivariateStatistics`"]

points = {{31.4799, 432.849}, {-195.826, 356.419}, {210.029, 
    121.779}, {76.1586, -4.47992}, {26.6061, -6.97711}, {160.236, 
    27.0398}, {203.248, 241.865}, {-225.351, 218.912}, {-159.899, 
    109.93}, {-106.465, 58.164}, {-1.98952*10^-11, 442.}, {-229.146, 
    324.489}, {-31.4799, 9.15114}, {195.826, 85.5807}, {-210.029, 
    320.221}, {-76.1586, 446.48}, {-26.6061, 448.977}, {-160.236, 
    414.96}, {-203.248, 200.135}, {225.351, 223.088}, {159.899, 
    332.07}, {106.465, 383.836}, {-1.98952*10^-11, 
    1.83076*10^-11}, {229.146, 117.511}};

Show[ListPlot[points, PlotMarkers -> {Automatic, Medium}], 
 Graphics[EllipsoidQuantile[points, {0.1, 1.}]], PlotRange -> All]

enter image description here

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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