Least square fit of ellipse worsens with increasing border thickness

Question

I use the least square method to fit ellipses to data points that originate from a binary image. The least square fit minimizes the sum of squared algebraic distances between data points and the closest point to a non-centered, rotated ellipse. For the definition of algebraic distance, see $[1]$ and following snippet:

In the article also $6$ lines of Matlab code are given that implements the least square method for an ellipse-only fit (never parabola or hyperbola are fitted "accidentially"). This code was used for the fit.

Below we see $4$ versions of the same ellipse where from left to right the border thickness was increased inwards and outwards by the same amount. The least square fit works well for thin borders but with increasing border thickness the fit systematically worsens. This is visible by the naked eye for the most right figure (data points for this ellipse can be downloaded here). The fitted ellipse is not in the center of the elliptical shaped data points. With increasing border thickness the large semi-axis decreases and the small semi-axis increases, the other fit parameter are in principle not affected.

The same fits I get if I use another algorithm suggested in $[2]$ (using the Matlab code included in the paper). Using Mathematica with another code delivers an even worse fit:

The inner and outer borders of the data points are approximately confined by ellipses (green, manually fitted). The ellipse that I would expect by a least square fit would lie centrally between the green ellipses and would be similiar to the manually drawn red ellipse.

How to explain and alleviate this phenomen to get good fits? Maybe there is a short algorithm/code? A good fit should go approximately through the center of the data points.

Further information

• there are no overlapping data points
• all data points have the same weight
• original data points lie on a grid and before fit Gaussian noise is added (standard deviation $\ll$ grid spacing)
• thinning thick borders in the original images is not the desired solution

References

$[1]$ Andrew Fitzgibbon, Maurizio Pilu, and Robert B. Fisher: Direct Least Square Fitting of Ellipses, Pattern Analysis and Machine Intelligence, 21, 476 (1999). PDF

$[2]$ Radim Halir and Jan Flusser: Numerically Stable Direct Least Squares Fitting of Ellipses, Proc. Sixth Int. Conf. in Central Europe on Computer Graphics and Visualization, 1, 125 (1998) PDF

Answer 1

You can make the least-squares method work, but you have to be careful about which least-squares system you solve.

Clearly, the equation $Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0$ for the ellipse isn’t unique: multiplying $A, B, C, D, E, F$ by the same constant gives you another equation for the same ellipse. So you can’t simply minimize

$$\sum (Ax^2 + Bxy + Cy^2 + Dx + Ey + F)^2$$

over $A, B, C, D, E, F$, because you’ll get $A = B = C = D = E = F = 0$; you need to add a normalizing constraint to exclude this solution. But if you add the wrong constraint, for example $F = 1$, you’ll bias the solution towards ellipses where $A, B, C, D, E$ are smaller relative to $F$.

The right constraint is $A + C = 1$, because $A + C$ is invariant over isometries of $(x, y)$. Minimize

$$\sum (Bxy + C(y^2 - x^2) + Dx + Ey + F + x^2)^2$$

over $B, C, D, E, F$, and then let $A = 1 - C$.

An advantage of this method over one using barycenters and inertial moments is that it still works well with a non-uniform distribution of points.

Python code for these figures:

import numpy as np
from matplotlib.patches import Ellipse
import matplotlib.pyplot as plt
with open("ellipse_data.txt") as file:
    points = np.array([[float(s) for s in line.split()] for line in file])
xs, ys = points.T
Parameters for Ax² + Bxy + Cy² + Dx + Ey + F = 0
B, C, D, E, F = np.linalg.lstsq(
    np.c_[xs * ys, ys ** 2 - xs ** 2, xs, ys, np.ones_like(xs)], -(xs ** 2)
)[0]
A = 1 - C
Parameters for ((x-x0)cos θ + (y-y0)sin θ)²/a² + (-(x-x0)sin θ + (y-y0)cos θ)²/b² = 1
M = np.array([[A, B / 2, D / 2], [B / 2, C, E / 2], [D / 2, E / 2, F]])
λ, v = np.linalg.eigh(M[:2, :2])
x0, y0 = np.linalg.solve(M[:2, :2], -M[:2, 2])
a, b = np.sqrt(-np.linalg.det(M) / (λ[0] * λ[1] * λ))
θ = np.arctan2(v[1, 0], v[0, 0])
ax = plt.axes(aspect="equal")
ax.add_patch(
    Ellipse((x0, y0), 2 * a, 2 * b, θ * 180 / np.pi, facecolor="none", edgecolor="red")
)
ax.scatter(xs, ys, s=0.5)
plt.show()

Answer 2

Using the symmetry, after a convenient rotation.

First we calculate the inertia matrix centered at the barycenter:

Given a set of points $S = \{s_k\} = \{x_k,y_k\},\ \ k=1,\cdots,n$ we calculate the barycenter $g = \frac 1n\sum_{k=1}^n s_k$. After that we calculate the inertia matrix

$$ I = \left( \begin{array}{cc} I_{xx} & -I_{xy} \\ -I_{xy} & I_{yy} \\ \end{array} \right) = \left( \begin{array}{cc} 510964. & 157771. \\ 157771. & 174852. \\ \end{array} \right) $$

with eigenvectors $v_1 = \{-0.929802,-0.36806\},\ \ v_2 = \{0.36806,-0.929802\}$

Then we rotate the centered data points $S-g$ by an angle given by $\alpha = -\arccos\left(v_1\cdot \vec e_y\right) = -1.94772$ giving $S_{\alpha}$

Now we define a distance to perform the fitting as:

$$ \delta(\rho_k,\theta_k,a,c) = \left|\rho_k - \frac{1}{\sqrt{\frac{\cos^2\theta}{a^2}+\frac{\sin^2\theta}{c^2}}}\right| $$

After transforming to polar form the rotated data $S_{\alpha}\to \{\rho_k,\theta_k\}$ we solve the minimization problem

$$ \min_{a,c}\sum_{k=1}^n \delta(\rho_k,\theta_k,a,c) $$

having the result

$$ E_{\alpha}=\frac{x^2}{33.959^2}+\frac{y^2}{12.0246^2}-1 = 0 $$

and recovering the original reference frame

$$ E=0.0060966 x^2+0.00414015 x y-0.393051 x+0.00168657 y^2-0.204216 y+6.60746 = 0 $$

NOTE

Attached a MATHEMATICA script to implement the process.

Clear["Global`*"]
dr[r_, t_] := r - 1/Sqrt[Cos[t]^2/a^2 + Sin[t]^2/c^2]
m0 = Import["/path_to_data/ellipse_data.txt", "Table"];
n = Length[m0];
g = Total[m0]/n;
{X, Y} = m0 // Transpose;
ixx = Sum[(Y[[k]] - g[[2]])^2, {k, 1, n}];
iyy = Sum[(X[[k]] - g[[1]])^2, {k, 1, n}];
ixy = -Sum[(X[[k]] - g[[1]]) (Y[[k]] - g[[2]]), {k, 1, n}];
mI = {{ixx, ixy}, {ixy, iyy}};
{Lambda, {v1, v2}} = Eigensystem[mI];
ang = -ArcCos[v1.{0, 1}/Norm[v1]];
m0r = Table[RotationMatrix[ang].(m0[[k]] - g), {k, 1, n}];
polardata = Table[{Norm[m0r[[k]]], ArcTan[m0r[[k, 1]], m0r[[k, 2]]]}, {k, 1, n}];
pobj = Sum[(dr[polardata[[k, 1]], polardata[[k, 2]]])^2, {k, 1, n}];
restrs = {0 < c < 15, 15 < a < 40};
psol = NMinimize[Join[{pobj}, restrs], {a, c}];
pellipse = x^2/a^2 + y^2/c^2 - 1 /. psol[[2]];
ellipser = pellipse /. {Thread[{x, y} -> RotationMatrix[ang].{XX - g[[1]], YY - g[[2]]}]} // N // Expand
grg = Graphics[{Red, PointSize[0.05], Point[g]}];
gr00 = ListPlot[m0, AspectRatio -> 2, PlotStyle -> {Thick, Black}];
gr3 = ContourPlot[ellipser == 0, {XX, 0, 45}, {YY, 0, 75},ContourStyle -> {Thick, Red}, PlotRange -> All];
Show[gr00, gr3, grg]

Answer 3

How is the objective function defined? In other words, what exactly is the quantity that is being minimized here? Just saying "least squares" does not adequately specify how the fit is performed.

For example, the fit might be performed on the following parameters:

1. The center $(x_0, y_0)$ of the ellipse
2. The semimajor and semiminor axis lengths $a, b$.
3. The rotation angle $\theta$ of the ellipse.

But then, given such an ellipse, how do we determine the error? Is it the sum of the squared distances of each data point to the closest point on the ellipse? Or is it the the sum of the squared vertical distances of each data point to the closest point on the ellipse? Before you question the latter option, note that this is precisely how a linear least-squares fit is performed: given observed $(x_i, y_i)$, we are finding $(a,b)$ that minimizes $$\sum_{i=1}^n (y_i - (ax_i + b))^2.$$