如何将2D椭圆拟合到给定点 [英] How to fit a 2D ellipse to given points

问题描述

我想通过椭圆函数拟合 2D数组:(x/a)²+(y/b)²= 1 ----> (因此得到a和b)

I would like to fit a 2D array by an elliptic function: (x / a)² + (y / b)² = 1 ----> (and so get the a and b)

然后，可以将其重新绘制在我的图表上. 我在互联网上发现了许多示例，但没有人知道这个简单的笛卡尔方程式.我可能搜索不好！ 我认为解决此问题的基本方法可以帮助很多人.

And then, be able to replot it on my graph. I found many examples on internet, but no one with this simple Cartesian equation. I probably have searched badly ! I think a basic solution for this problem could help many people.

以下是数据示例:

遗憾的是，我无法输入值...因此，假设我有一个X，Y数组，用于定义每个点的坐标.

Sadly, I can not put the values... So let's assume that I have an X,Y arrays defining the coordinates of each of those points.

推荐答案

可以使用最小二乘法直接解决.您可以将其构架为使数量的平方和最小化(alpha * x_i ^ 2 + beta * y_i ^ 2-1-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是所有列的列向量.

This can be solved directly using least squares. You can frame this as minimizing the sum of squares of quantity (alpha * x_i^2 + beta * y_i^2 - 1) where alpha is 1/a^2 and beta is 1/b^2. You have all the x_i's in X and the y_i's in Y so you can find the minimizer of ||Ax - b||^2 where A is an Nx2 matrix (i.e. [X^2, Y^2]), x is the column vector [alpha; beta] and b is column vector of all ones.

下面的代码解决了形式为Ax ^ 2 + Bxy + Cy ^ 2 + Dx + Ey = 1的椭圆的更普遍的问题，尽管这个想法是完全相同的.打印语句给出0.0776x ^ 2 + 0.0315xy + 0.125y ^ 2 + 0.00457x + 0.00314y = 1并且生成的椭圆图像也在下面

The following code solves the more general problem for an ellipse of the form Ax^2 + Bxy + Cy^2 + Dx +Ey = 1 though the idea is exactly the same. The print statement gives 0.0776x^2 + 0.0315xy+0.125y^2+0.00457x+0.00314y = 1 and the image of the ellipse generated is also below

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椭圆拟合到给定点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！