Как подогнать 2D-эллипс к заданным точкам

Это можно решить напрямую с помощью метода наименьших квадратов. Вы можете сформулировать это как минимизацию суммы квадратов количества (альфа * x_i^2 + бета * y_i^2 - 1), где альфа равна 1/a^2, а бета равна 1/b^2. У вас есть все x_i в X и y_i в Y, поэтому вы можете найти минимизатор ||Ax - b||^2, где A — матрица Nx2 (т. е. [X^2, Y^2]), x — это вектор-столбец [альфа; beta], а b — вектор-столбец всех единиц.

Следующий код решает более общую задачу для эллипса формы Ax^2 + Bxy + Cy^2 + Dx +Ey = 1, хотя идея точно такая же. Оператор печати дает 0,0776x ^ 2 + 0,0315xy + 0,125y ^ 2 + 0,00457x + 0,00314y = 1, а изображение сгенерированного эллипса также ниже.

import numpy as np
import matplotlib.pyplot as plt
alpha = 5
beta = 3
N = 500
DIM = 2

np.random.seed(2)

# Generate random points on the unit circle by sampling uniform angles
theta = np.random.uniform(0, 2*np.pi, (N,1))
eps_noise = 0.2 * np.random.normal(size=[N,1])
circle = np.hstack([np.cos(theta), np.sin(theta)])

# Stretch and rotate circle to an ellipse with random linear tranformation
B = np.random.randint(-3, 3, (DIM, DIM))
noisy_ellipse = circle.dot(B) + eps_noise

# Extract x coords and y coords of the ellipse as column vectors
X = noisy_ellipse[:,0:1]
Y = noisy_ellipse[:,1:]

# Formulate and solve the least squares problem ||Ax - b ||^2
A = np.hstack([X**2, X * Y, Y**2, X, Y])
b = np.ones_like(X)
x = np.linalg.lstsq(A, b)[0].squeeze()

# Print the equation of the ellipse in standard form
print('The ellipse is given by {0:.3}x^2 + {1:.3}xy+{2:.3}y^2+{3:.3}x+{4:.3}y = 1'.format(x[0], x[1],x[2],x[3],x[4]))

# Plot the noisy data
plt.scatter(X, Y, label='Data Points')

# Plot the original ellipse from which the data was generated
phi = np.linspace(0, 2*np.pi, 1000).reshape((1000,1))
c = np.hstack([np.cos(phi), np.sin(phi)])
ground_truth_ellipse = c.dot(B)
plt.plot(ground_truth_ellipse[:,0], ground_truth_ellipse[:,1], 'k--', label='Generating Ellipse')

# Plot the least squares ellipse
x_coord = np.linspace(-5,5,300)
y_coord = np.linspace(-5,5,300)
X_coord, Y_coord = np.meshgrid(x_coord, y_coord)
Z_coord = x[0] * X_coord ** 2 + x[1] * X_coord * Y_coord + x[2] * Y_coord**2 + x[3] * X_coord + x[4] * Y_coord
plt.contour(X_coord, Y_coord, Z_coord, levels=[1], colors=('r'), linewidths=2)

plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

Следуя предложению ErroriSalvo, вот полный процесс подбора эллипса с помощью SVD. Массивы x, y — это координаты заданных точек, допустим, точек N. Затем U, S, V получаются из СВД центрированного координатного массива формы (2, N). Итак, U — ортогональная матрица 2 на 2 (поворот), S — вектор длины 2 (сингулярные значения), а V, который нам не нужен, — ортогональная матрица N на N.

Линейная карта, преобразующая единичный круг в эллипс наилучшего соответствия, представляет собой

sqrt(2/N) * U * diag(S) 

где diag(S) — диагональная матрица с сингулярными значениями на диагонали. Чтобы понять, зачем нужен множитель sqrt(2/N), представьте, что точки x, y взяты равномерно из единичного круга. Тогда sum(x**2) + sum(y**2) равно N, поэтому матрица координат состоит из двух ортогональных строк длины sqrt(N/2), следовательно, ее норма (наибольшая сингулярная величина) равна sqrt(N/2). Нам нужно уменьшить это значение до 1, чтобы получить единичный круг.

N = 300
t = np.linspace(0, 2*np.pi, N)
x = 5*np.cos(t) + 0.2*np.random.normal(size=N) + 1
y = 4*np.sin(t+0.5) + 0.2*np.random.normal(size=N)
plt.plot(x, y, '.')     # given points

xmean, ymean = x.mean(), y.mean()
x -= xmean
y -= ymean
U, S, V = np.linalg.svd(np.stack((x, y)))

tt = np.linspace(0, 2*np.pi, 1000)
circle = np.stack((np.cos(tt), np.sin(tt)))    # unit circle
transform = np.sqrt(2/N) * U.dot(np.diag(S))   # transformation matrix
fit = transform.dot(circle) + np.array([[xmean], [ymean]])
plt.plot(fit[0, :], fit[1, :], 'r')
plt.show()

Вы можете попробовать разложение матрицы данных по сингулярным значениям. https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.svd.html

Сначала центрируйте данные, вычитая средние значения X, Y из каждого столбца соответственно.

X=X-np.mean(X)
Y=Y-np.mean(Y)

D=np.vstack(X,Y)

Затем примените SVD и извлеките -собственные значения (члены s) -> длина оси -собственные векторы (U) -> ориентация оси

U, s, V = np.linalg.svd(D, full_matrices=True)

Это должно быть методом наименьших квадратов. Конечно, все может быть сложнее, см. https://www.emis.de/journals/BBMS/Bulletin/sup962/gander.pdf