Cómo ajustar una elipse a datos 2D

Me gustaría encontrar el "mejor" ajuste de una elipse a formas contiguas, posiblemente cóncavas como:

Que he probado

Pensé que se podía asignar la dirección de los ejes mayor y menor. $\vec a, \vec b$de la elipse mapeando los valores de píxel a las coordenadas, restando la media y guardando los dos vectores propios más grandes de un PCA. Parece que funciona bastante bien para encontrar la dirección:

Mi problema es determinar la longitud de estos dos vectores. Por ahora, he usado$\sqrt \lambda_1$ $\sqrt \lambda_2$, de los valores propios de la PCA. Esto parece subestimar la longitud. ¿Cómo puedo determinar$|\vec a|$y $|\vec b|$ o como alternativa, ¿se ajusta mejor una elipse a estas formas?

Continuando con mi respuesta eliminada ... Si toma una elipse llena y proyecta todos los puntos en el $x$eje, se proyectarán más puntos cerca del origen que en los extremos, en una distribución de forma circular. No es una distribución gaussiana, y no la distribución uniforme que mencioné en la analogía 1-D en mi respuesta eliminada. La distribución resultante en realidad tiene pdf$p(x) = \sqrt{(1 - (\frac{x}{r})^2)}$, y desde allí puede calcular que la desviación estándar es $\frac{r}{2}$.

Por lo tanto, si los datos se distribuyen uniformemente en el interior de una elipse de radios $a, b$ (cuyos ejes son el $x$ y $y$ ejes), la desviación estándar de la $x$ coordenada es $\frac{a}{2}$ y de la $y$ coordenada es $\frac{b}{2}$. Entonces, el factor de corrección que necesita usar es simplemente 2.

Aquí hay un ejemplo trabajado en python que recupera el centro (matriz de traducción), la matriz de rotación y los radios de una elipse desde puntos muestreados aleatoriamente desde su 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)

Compruebe el ajuste de mínimos cuadrados directos de elipses por Fitzgibbon et al. ¡Es un problema de valor propio simple, cuyo tamaño no es proporcional al número de píxeles en su muestra! El único paso que depende de la cantidad de píxeles que le arroje es el cálculo de una matriz de dispersión, que aún es$O(n)$.