如何将 2D 椭圆拟合到给定点

Question

我想用椭圆函数拟合 2D 数组：(x / a)² + (y / b)² = 1 ----> (所以得到a和b)

然后，能够在我的图表上重新绘制它。我在互联网上找到了很多例子，但没有人知道这个简单的笛卡尔方程。我可能搜索得很糟糕！ 我认为这个问题的基本解决方案可以帮助很多人。

以下是数据示例：

遗憾的是，我无法输入值...所以让我们假设 我有一个 X、Y 数组定义每个点的坐标。

Answer 1

您可以尝试对数据矩阵进行奇异值分解。 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

Answer 2

根据ErroriSalvo的建议，下面是使用SVD拟合椭圆的完整过程。数组 x, y 是给定点的坐标，假设有 N 个点。然后U，S，V由形状为（2，N）的中心坐标阵列的SVD得到。因此，U 是一个 2 x 2 正交矩阵（旋转），S 是一个长度为 2 的向量（奇异值），而我们不需要的 V 是一个 N x 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()

但是如果你假设没有旋转，那么np.sqrt(2/N) * S就是你所需要的；这些是椭圆方程中的 a 和 b。

Answer 3

这可以直接用最小二乘法求解。您可以将其定义为最小化数量平方和 (alpha * x_i^2 + beta * y_i^2 - 1)，其中 alpha 为 1/a^2，beta 为 1/b^ 2. X 中有所有 x_i，Y 中有 y_i，因此您可以找到 ||Ax - b||^2 的最小值，其中 A 是 Nx2 矩阵（即 [X ^2, Y^2]), x为列向量[alpha; 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()

如何将 2D 椭圆拟合到给定点

How to fit a 2D ellipse to given points

python-2.7

numpy

cartesian