将椭圆拟合到一组二维点

Question

我正在尝试将椭圆拟合到由 x 和 y 坐标描述的一组点。

我在这里找到了详细的说明http://nicky.vanforeest.com/misc/fitEllipse/fitEllipse.html 并尝试了代码，但它似乎不起作用。

它正确地找到了中心，但是角度和轴完全错误，如下图所示： https://i.imgur.com/VLEeNKQ.png

红色点是我的数据点，蓝色点是从获得的参数中绘制的椭圆。现在，数据不是一个完美的椭圆，但拟合度远不及它。我想让拟合更接近实际数据。

这是有问题的代码。

import numpy as np
from numpy.linalg import eig, inv

def fitEllipse(x,y):
    x = x[:,np.newaxis]
    y = y[:,np.newaxis]
    D =  np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
    S = np.dot(D.T,D)
    C = np.zeros([6,6])
    C[0,2] = C[2,0] = 2; C[1,1] = -1
    E, V =  eig(np.dot(inv(S), C))
    n = np.argmax(np.abs(E))
    a = V[:,n]
    return a

def ellipse_center(a):
    b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
    num = b*b-a*c
    x0=(c*d-b*f)/num
    y0=(a*f-b*d)/num
    return np.array([x0,y0])


def ellipse_angle_of_rotation( a ):
    b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
    return 0.5*np.arctan(2*b/(a-c))


def ellipse_axis_length( a ):
    b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
    up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
    down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
    down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
    res1=np.sqrt(up/down1)
    res2=np.sqrt(up/down2)
    return np.array([res1, res2])

def ellipse_angle_of_rotation2( a ):
    b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
    if b == 0:
        if a > c:
            return 0
        else:
            return np.pi/2
    else: 
        if a > c:
            return np.arctan(2*b/(a-c))/2
        else:
            return np.pi/2 + np.arctan(2*b/(a-c))/2

这是我的完整数据集。有谁知道为什么它没有正确安装？

# --------------------------------------------------------------------------
x = np.array([ 5727.53135,  7147.62235, 10330.93573,  8711.17228, 7630.40262,
        4777.24983,  4828.27655,  9449.94416,  5203.81323,  6299.44811,
        6494.21906])

y = np.array([67157.77567 , 66568.50068 , 55922.56257 , 54887.47348 ,
       65150.14064 , 66529.91705 , 65934.25548 , 55351.57612 ,
       63123.5103  , 67181.141725, 56321.36025 ])
# -----------------------------------------------------------------------------


a = fitEllipse(x,y)
center = ellipse_center(a)
#phi = ellipse_angle_of_rotation(a)
phi = ellipse_angle_of_rotation2(a)
axes = ellipse_axis_length(a)

# get the individual axes
a, b = axes

from matplotlib.patches import Ellipse
import matplotlib.pyplot as plt

ell = Ellipse(center, a, b, phi)

fig, ax = plt.subplots(subplot_kw={'aspect': 'equal'})
ax.add_artist(ell)
ell.set_clip_box(ax.bbox)
ax.set_xlim(0, 100000)
ax.set_ylim(0, 100000)
plt.show()
scat = plt.scatter(x, y, c = "r")

Answer 1

您的代码绝对没问题，正是 patch 定义在这里造成了麻烦。 Ellipse 的 a 和 b 是全宽。所以你必须将你的拟合结果乘以2。此外，角度是度数，所以你必须乘以180/np.pi。最后，零不在同一个位置，所以要加上90。

长话短说

Ellipse(center, a, b, phi)

到

ell = Ellipse(center, 2 * a, 2 * b,  phi * 180 / np.pi + 90 ) 

你很好。