3
$\begingroup$

Given a set of $n$ two-dimensional points $\{(x_i, y_i), i=1,.., n\}$ I would like to find the ellipse that best fits them. The easiest approach I found was what is called the "Equation Error Model", in which we are trying to fit the points to the conic section model:

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

The coefficients $A,B,C,D,E,F$ are not unique because they can be scaled by any real non-zero number. Also, for an ellipse, $A$ cannot be zero. Therefore, we can choose to fix $A = 1$, then we would have to determine the other $5$ coefficients.

Accumulating the data from every point, we will have the linear system

$ G a = b $

where $ a = [B, C, D, E, F]^T$, and

$ G = \begin{bmatrix} x_1 y_1 && y_1^2 && x_1 && y_1 && 1 \\ x_2 y_2 && y_2^2 && x_2 && y_2 && 1 \\ x_3 y_3 && y_3^2 && x_3 && y_3 && 1 \\ \vdots \\ x_n y_n && y_n^2 && x_n && y_n && 1 \end{bmatrix}\hspace{5pt}, \hspace{25pt} b = \begin{bmatrix} - x_1^2 \\ - x_2^2 \\ - x_3^2 \\ \vdots \\ - x_n^2 \end{bmatrix} $

To solve this linear system for $a$, pre-multiply by $G^T$, you get

$ G^T G a = G^T b $

and finally,

$ a = (G^T G)^{-1} G^T b $

I simulated the above method, and it gave excellent results. My question is whether there are alternative methods to approach this best fitting ellipse problem. Your hints, suggestions, and solutions are much appreciated.

$\endgroup$
9
  • $\begingroup$ By using the pseudo inverse as you have, you are minimizing an L2 norm. This works great if your nose is Gaussian. If you have spurious points, most of the noise is reasonable but few points are not, then you can use RANSAC to address that. $\endgroup$ Commented May 28, 2022 at 14:48
  • 1
    $\begingroup$ Thank you. I'll look into the RANSAC algorithm. $\endgroup$ Commented May 28, 2022 at 15:45
  • $\begingroup$ scipython.com/blog/… $\endgroup$ Commented May 29, 2022 at 1:50
  • $\begingroup$ If that's the case, why bother editing it ? $\endgroup$ Commented Sep 9, 2022 at 9:41
  • $\begingroup$ This is not written without effort, my friend. $\endgroup$ Commented Sep 9, 2022 at 9:44

2 Answers 2

Reset to default
4
+50
$\begingroup$

It’s better to fix $A + C = 1$ rather than $A = 1$, for reasons I explained in Least square fit of ellipse worsens with increasing border thickness: $A + C$ is invariant under isometries of $(x, y)$, while $A$ is not.

Here’s a diagram using the same data from that question, comparing the ellipse correctly computed with $A + C = 1$ in red vs. incorrectly computed with $A = 1$ in orange. If you rotated the points, you’d find that the red ellipse matches their rotation exactly, while the orange ellipse wobbles around a bit.

diagram

$\endgroup$
1
$\begingroup$

The solution you state using the Moore-Penrose pseudoinverse coincides with the solution of the convex optimization problem:

$$\min_{a\in\mathbb{R}^5}||Ga-b||_2.$$

However, it is possible to use different norms for finding the "best" vector $a$. You can use, for instance, the $1$-norm which leads to the following convex optimization problem

$$\min_{a\in\mathbb{R}^5}||Ga-b||_1.$$

Note, however, that this problem does not necessarily have a closed-form optimal solution as in the $2$-norm case and one has to rely to numerical optimization methods for finding the minimizer. Using the $1$-norm has certain advantages over the $2$-norm such as "ignoring outliers" in a better way as it promotes sparsity, that is, it will try to set as many entries of $Ga-b$ as possible to zero.

You may also try the $\infty$-norm as well for fun and solve

$$\min_{a\in\mathbb{R}^5}||Ga-b||_\infty.$$

Those optimizations problems are all convex. In particular, minimizing the $1$-norm and the $\infty$-norm are linear programming problems. They can be solved using standard libraries in Python or Matlab, for instance.

$\endgroup$

You must log in to answer this question.