椭圆拟合以确定旋转（Python）

Question

我希望将椭圆拟合到我拥有的一些数据点。

我想要什么：使用椭圆来确定我的数据的旋转角度

数据：我拥有的数据是极坐标（θ，r）

theta = [0.0, 0.103, 0.206, 0.309, 0.412, 0.515, 0.618, 0.721, 0.824, 0.927, 1.03, 1.133, 1.236, 1.339, 1.442, 1.545, 1.648, 1.751, 1.854, 1.957, 2.06, 2.163, 2.266, 2.369, 2.472, 2.575, 2.678, 2.781, 2.884, 2.987, 3.09, 3.193, 3.296, 3.399, 3.502, 3.605, 3.708, 3.811, 3.914, 4.017, 4.12, 4.223, 4.326, 4.429, 4.532, 4.635, 4.738, 4.841, 4.944, 5.047, 5.15, 5.253, 5.356, 5.459, 5.562, 5.665, 5.768, 5.871, 5.974, 6.077, 6.18]

r = [84.48, 83.11, 77.67, 76.62, 90.12, 89.64, 84.07, 95.21, 104.63, 119.19, 125.19, 140.25, 146.33, 145.11, 164.0, 202.87, 214.81, 258.5, 281.94, 268.5, 224.76, 238.61, 270.08, 245.86, 220.04, 179.98, 181.51, 189.53, 172.87, 153.29, 138.32, 156.67, 146.21, 129.28, 139.76, 132.12, 138.73, 133.83, 136.15, 172.02, 163.2, 157.6, 142.73, 130.79, 130.24, 128.88, 124.7, 119.37, 115.28, 118.02, 117.89, 121.73, 115.13, 103.02, 84.43, 83.69, 82.26, 87.87, 88.84, 92.53, 94.67]

目前的算法：

1. 定义残差和残差的雅可比
2. 使用 scipy optimize.leastsq

（这里是有兴趣的人的演练https://scipython.com/book/chapter-8-scipy/examples/non-linear-fitting-to-an-ellipse/）

然而，在我的数据集上，偏心度是负数，如果它是一个椭圆 (0

我尝试添加一个取决于 theta 但到目前为止没有任何运气的轮换项。这是适合椭圆的代码，没有我的额外术语会搞砸一切：

import numpy as np
from scipy import optimize
import pylab

def f(theta, p):
    a, e = p
    return a * (1 - e**2)/(1 - e*np.cos(theta))

def residuals(p, r, theta):
    """ Return the observed - calculated residuals using f(theta, p). """
    return r - f(theta, p)

def jac(p, r, theta):
    """ Calculate and return the Jacobian of residuals. """
    a, e = p
    da = (1 - e**2)/(1 - e*np.cos(theta))
    de = (-2*a*e*(1-e*np.cos(theta)) + a*(1-e**2)*np.cos(theta))/(1 -
                                                        e*np.cos(theta))**2
    return -da,  -de
    return np.array((-da, -de)).T

def fit_ellipse(theta, r, p0 = (1,0.5)):
    # Initial guesses for a, e
    p0 = (1, 0.5)
    plsq = optimize.leastsq(residuals, p0, Dfun=jac, args=(r, theta), col_deriv=True)
    #return plsq
    print(plsq)
    
    pylab.polar(theta, r, 'o')
    theta_grid = np.linspace(0, 2*np.pi, 200)
    pylab.polar(theta_grid, f(theta_grid, plsq[0]), lw=2)
    pylab.show()

fit_ellipse(theta, r, p0 = (1,0.5))

Answer 1

不需要非线性回归（如果不需要特定的拟合标准）。一个简单的线性回归导致以下结果：

符号和符号与方程式一致。 (15-23) 来自https://mathworld.wolfram.com/Ellipse.ht

另外：回复评论。

用于绘制椭圆图形的方程是：

另一种绘制椭圆的方式（避免复杂的根）是：

参数 theta 必须从 0 到 2pi。