How to fit an ellipse to 2D data

Question

I'd like to find the "best" fit of an ellipse to contiguous, possibly concave shapes such as:

What have I tried?

I thought that one could assign the direction of the major and minor axes $\vec a, \vec b$ of the ellipse by mapping the pixel values to coordinates, mean subtracting, and saving the largest two eigenvectors from a PCA. The seems to work rather well in finding the direction:

My problem is determining the length of these two vectors. For now, I've used $\sqrt \lambda_1$ $\sqrt \lambda_2$, from the eigenvalues of the PCA. This seems to underestimate the length. How can I determine $|\vec a|$, and $|\vec b|$ or alternatively, best fit an ellipse to these shapes?

Answer 1

Following up on my deleted answer... If you take a filled ellipse and project all the points onto the $x$ axis, more points will be projected near the origin than on the extremes, in a circular-shaped distribution. Not a gaussian distribution, and not the uniform distribution I mentioned in the 1-D analogy in my deleted answer. The resulting distribution actually has pdf $p(x) = \sqrt{(1 - (\frac{x}{r})^2)}$, and from there you can compute that the standard deviation is $\frac{r}{2}$.

Thus, if data is uniformly distributed in the interior of an ellipse of radii $a, b$ (whose axes are the $x$ and $y$ axes), the standard deviation of the $x$ coordinate is $\frac{a}{2}$ and of the $y$ coordinate is $\frac{b}{2}$. So the correction factor you need to use is simply 2.

Here is a worked example in python recovering the center (translation matrix), rotation matrix, and radii of an ellipse from points randomly sampled from its interior:

import numpy

# Generate some points distributed uniformely along an ellipse
x = 2 * (numpy.random.rand(10000, 1) - 0.5)
y = 2 * (numpy.random.rand(10000, 1) - 0.5)
d = (x / 0.5) ** 2 + (y / 0.25) ** 2
inside = numpy.where(d < 1)[0]
x = x[inside]
y = y[inside]
data = numpy.hstack((x, y)).T

# Now rotate by 0.5 rad and translate it to (4, -8)
angle = 0.5
rotation = numpy.array([
    [numpy.cos(0.4), -numpy.sin(0.4)],
    [numpy.sin(0.4), numpy.cos(0.4)]])

data = numpy.dot(rotation, data)
data[0, :] += 4
data[1, :] -= 8

# Step 1: center the data to get the translation vector.
print 'Translation', data.mean(axis=1)
data -= numpy.reshape(data.mean(axis=1), (2, 1))

# Step 2: perform PCA to find rotation matrix.
scatter = numpy.dot(data, data.T)
eigenvalues, transform = numpy.linalg.eig(scatter)
print 'Rotation matrix', transform

# Step 3: Rotate back the data and compute radii.
# You can also get the radii from smaller to bigger
# with 2 / numpy.sqrt(eigenvalues)
rotated = numpy.dot(numpy.linalg.inv(transform), data)
print 'Radii', 2 * rotated.std(axis=1)

Answer 2

Check Direct Least Squares Fitting of Ellipses by Fitzgibbon et al. It's a simple eigenvalue problem, the size of which is not proportional to the number of pixels in your sample! The only step dependent on the number of pixels you throw at it is the computation of a scatter matrix, which is still $O(n)$.