Least-squares ellipse fitting

Question

I am trying to find a least-squares ellipse fit for a set of 100 data points $(x,y)$.

Now I have found the values of $A,B,C,D,E,F$ according to the conical equation of the ellipse $$ Ax^2+Bxy+Cy^2+Dx+Ey+F=0 $$ I would like to know how to find the points that actually lie on this ellipse. From my basic understanding, if I substitute a value of $x$ in the above equation, it should give me the corresponding value of $y$.

When I do the above, I get a straight line and not really a fitted ellipse. How can I find the fitted ellipse?

My task is to plot these points so that I can see the best possible fit. For reference see [link]. This is the source of ellipse fitting that I am currently using.

I appreciate help from anyone who has experience with this. I am sorry if I am lacking some basic mathematical knowledge, but from what I understand, it isn't all that straightforward.

Regards

Arj

Answer 1

Given a ellipse as

$$ E(x,y)=Ax^2+B xy+Cy^2+D x+ E y + F=0 $$

and a data set as $(x_k,y_k),\{k=1\cdots,n\}$ a typical fitting process is to minimize the residuals at each point, squared so defining

$$ R(A,B,C,D,E,F) = \sum_{k=1}^nE(x_k,y_k)^2 $$

the problem reads as

$$ \min_{A,B,C,D,E,F}R\ \ \ \ \text{s. t.}\ \ \ \ \cases{B^2\lt 4A C\\ A^2+B^2\gt \mu} $$

those restrictions are needed first to guarantee that the fitted conic is an ellipse and second to avoid the undesirable trivial minimum $A=B=C=D=E=F=0$. Follows a MATHEMATICA script showing the procedure

(** DATA GENERATION **)
Clear[a, b, c, d, e, f, ellipse]
r[a_, b_, theta_, t_] := a Sqrt[2]/Sqrt[(1 + (a/b)^2) + (1 - (a/b)^2) Cos[2 (t - theta)]]; 
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]}]];
gr1 = ListPlot[data, PlotStyle -> {Red,PointSize[Medium]}];
(** DATA FITTING **)
M = {{aa, bb}, {bb, cc}};
B = {dd, ee};
ellipse[x_, y_] := {x, y}.M.{x, y} + B.{x, y} + ff
res[k_] := ellipse[data[[k, 1]], data[[k,2]]]^2
sol = Minimize[{Sum[res[k], {k, 1,Length[data]}], 4 aa cc > bb^2, aa^2 + bb^2 > 0.3}, {aa, bb, cc, dd, ee, ff}]
ellipse0 = ellipse[x, y] /. sol[[2]]
gr2 = ContourPlot[ellipse0 == 0, {x, 0, 2},{y, 0, 2}, ContourStyle -> {Thick, Blue}];
Show[gr1, gr2, AspectRatio -> 1]

Answer 2

A sure way to determine whether a point lies on the ellipse is to substitute the point's x- and y-coordinates into the equation and see whether the equation is exactly satisfied. (Your method is also valid, but requires more work.)

An example of how to use MATLAB to plot a curve and points is at https://www.mathworks.com/help/curvefit/fit.html. You will need to solve your equation for y. Be careful, though, because a point may appear to be on the ellipse but actually not be on it.

Answer 3

Define the error associated with selecting a particular set of $A,B,C,D,E,F$ to be

$ \varepsilon_i = A x_i^2 + B x_i y_i + C y_i^2 + D x_i + E y_i + F $

This is the error for the $i$-th point $(x_i, y_i)$.

We want to minimize

$ J = \displaystyle \sum_{i=1}^N \varepsilon_i^2 $

To make this work we need to somehow normalize the set $A,B,C,D,E,F$. A good way to do this is to impose that

$ A + C = 1$

i.e. we take $C = 1 - A $

The error function is now

$J = \displaystyle \sum_{i=1}^N \left( A x_i^2 + B x_i y_i + (1 - A) y_i^2 + D x_i + E y_i + F \right)^2 $

At the minimum of $E$, its partial derivatives with respect to $A, B, D, E, F$ are zero.

This leads to the following "normal" equations:

$ \dfrac{\partial J}{\partial A} = \displaystyle \sum_{i=1}^N 2 \left(A x_i^2 + B x_i y_i + (1 - A) y_i^2 + D x_i + E y_i + F \right) ( x_i^2 - y_i^2 ) = 0 $

$ \dfrac{\partial J}{\partial B} = \displaystyle \sum_{i=1}^N 2 \left(A x_i^2 + B x_i y_i + (1 - A) y_i^2 + D x_i + E y_i + F \right) ( x_i y_i ) = 0 $

$ \dfrac{\partial J}{\partial D} = \displaystyle \sum_{i=1}^N 2 \left(A x_i^2 + B x_i y_i + (1 - A) y_i^2 + D x_i + E y_i + F \right) ( x_i ) = 0 $

$ \dfrac{\partial J}{\partial E} = \displaystyle \sum_{i=1}^N 2 \left(A x_i^2 + B x_i y_i + (1 - A) y_i^2 + D x_i + E y_i + F \right) ( y_i ) = 0 $

$ \dfrac{\partial J}{\partial F} = \displaystyle \sum_{i=1}^N 2 \left(A x_i^2 + B x_i y_i + (1 - A) y_i^2 + D x_i + E y_i + F \right) = 0 $

The above are $5$ non-homogeneous linear equations in $A, B, D, E, F$, that can be solved readily for the $5$ unknowns.

We're done.

To obtain the standard parameters of this ellipse, so that you can plot it, follow the steps below.

First, note that

$ A x^2 + B x y + C y^2 + D x + E y + F = 0 $

Can be written as

$r^T Q r + b^T r + F = 0 $

where

$ r = \begin{bmatrix} x \\ y \end{bmatrix} $

$ Q = \begin{bmatrix} A && \dfrac{B}{2} \\ \dfrac{B}{2} && C \end{bmatrix} $

$ b = \begin{bmatrix} D \\ E \end{bmatrix} $

The next step is to find the center of the ellipse, it is given by

$ r_0 = -\dfrac{1}{2} Q^{-1} b = -\dfrac{1}{2} \begin{bmatrix} A && \dfrac{B}{2} \\ \dfrac{B}{2} && C \end{bmatrix}^{-1} \begin{bmatrix} D \\ E \end{bmatrix} $

With this, the ellipse equation becomes

$ (r - r_0)^T Q (r - r_0) = - F + r_0^T Q r_0 $

Define the constant $K = - F + r_0^T Q r_0 $. Divide both sides of the equation by this constant. This will give you

$ (r - r_0)^T Q' (r - r_0) = 1 $

which is the standard algebraic equation of an ellipse.

The next step is to diagonalize $Q'$ by finding its two eigenvalues and the corresponding unit eigenvectors. Then you can write $Q'$ as follows

$ Q' = R D R^T $

where $D $ is the diagonal matrix whose diagonal entries are the two eigenvalues, and matrix $R$ has the two corresponding unit eigenvectors as its columns.

Equation of the ellipse becomes

$ (r - r_0)^T R D R^T (r - r_0) = 1 $

Introduce the new vector $p = D^{(1/2)} R^T (r - r_0) $

Where $D^{(1/2)} $ is the square root of $D$, i.e., it is the diagonal matrix whose diagonal entries are the square roots of the diagonal entries of $D$.

It follows that $r = r_0 + R D^{-(1/2)} p $

Plugging this into the above equation, we get

$ p^T p = 1 $

which means that $ p $ is a unit vector. The standard unit vector in $\mathbb{R}^2 $ is

$ p(t) = \begin{bmatrix} \cos t \\ \sin t \end{bmatrix} $

Therefore, finally, we now have an explicit parametrization of the ellipse, which is

$ r = r_0 + R D^{-(1/2)} p(t) $