将椭圆拟合到python中的一组数据点

Posted 2023-03-12

【中文标题】将椭圆拟合到python中的一组数据点

【英文标题】：Fitting an ellipse to a set of data points in python

【发布时间】：2017-02-03 06:25:40

【问题描述】：

我有一个二维点 (x,y)，我想使用这篇文章来拟合椭圆

fit a ellipse in Python given a set of points xi=(xi,yi)

但我的结果是axes = [ 0.93209407 nan]，因为在函数ellipse_axis_length 中down2 是一个负数，所以res2 无效，怎么办？如果我想根据数据集画一个椭圆，并计算数据点和椭圆之间的误差，我该怎么办？

代码是这样的：

import numpy as np
import numpy.linalg as linalg
import matplotlib.pyplot as plt

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 =  linalg.eig(np.dot(linalg.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 find_ellipse(x, y):
    xmean = x.mean()
    ymean = y.mean()
    x = x - xmean
    y = y - ymean
    a = fitEllipse(x,y)
    center = ellipse_center(a)
    center[0] += xmean
    center[1] += ymean
    phi = ellipse_angle_of_rotation(a)
    axes = ellipse_axis_length(a)
    x += xmean
    y += ymean
    return center, phi, axes

if __name__ == '__main__':

    points = [( 0 , 3),
        ( 1 , 2),
        ( 1 , 7),
        ( 2 , 2),
        ( 2 , 4),
        ( 2 , 5),
        ( 2 , 6),
        ( 2 ,14),
        ( 3 , 4),
        ( 4 , 4),
        ( 5 , 5),
        ( 5 ,14),
        ( 6 , 4),
        ( 7 , 3),
        ( 7 , 7),
        ( 8 ,10),
        ( 9 , 1),
        ( 9 , 8),
        ( 9 , 9),
        (10,  1),
        (10,  2),
        (10 ,12),
        (11 , 0),
        (11 , 7),
        (12 , 7),
        (12 ,11),
        (12 ,12),
        (13 , 6),
        (13 , 8),
        (13 ,12),
        (14 , 4),
        (14 , 5),
        (14 ,10),
        (14 ,13)]

    fig, axs = plt.subplots(2, 1, sharex = True, sharey = True)
    a_points = np.array(points)
    x = a_points[:, 0]
    y = a_points[:, 1]
    axs[0].plot(x,y)
    center, phi, axes = find_ellipse(x, y)
    print "center = ",  center
    print "angle of rotation = ",  phi
    print "axes = ", axes

    axs[1].plot(x, y)
    axs[1].scatter(center[0],center[1], color = 'red', s = 100)
    axs[1].set_xlim(x.min(), x.max())
    axs[1].set_ylim(y.min(), y.max())

    plt.show()

【问题讨论】：

看看the original post可能会有帮助

我认为您的数据量可能太小（

【参考方案1】：

我认为代码有错误。

我发现了一些其他问题（1、2）关于将椭圆拟合到一组数据点，它们都使用来自here 的同一段代码。

在装修中：

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

使用E 中的最大绝对特征值选择特征值-特征向量对。然而，在 Fitzgibbon、Pilu 和 Fischer 在 Fitzgibbon, A.W.、Pilu, M. 和 Fischer R.B. 的原始论文中，椭圆的直接最小二乘拟合，1996：

我们注意到特征系统 (6) 的解给出了 6 特征值-特征向量对 (\lambda_i, u_i)。这些对中的每一对 如果 (8) 中的平方根下的项，则产生局部最小值 是积极的。一般来说，S 是正定的，所以分母 u_i^T S u_i 对所有 u_i 都是积极的。因此平方根 如果\lambda_i>0 则存在。

他们进一步证明了只有一个正的特征值解，这也是6个特征值中的最大值。

但是，通过np.argmax(np.abs(E)) 找到这个最大值可以选择一个原本为负的值，从而给出错误的答案。

我找到了一个可以说明问题的具体示例。下面是一个 (x,y) 坐标数组：

159.598484728,45.0095214844
157.713012695,45.7333132048
156.163772773,46.6618041992
154.15536499,47.3460795985
152.140428382,48.0673522949
150.045213426,48.4620666504
148.194464489,47.868850708
145.55770874,47.6356541717
144.0753479,48.6449446276
144.19475866,50.200668335
144.668289185,51.7677851197
145.55770874,53.033584871
147.632995605,53.5380252111
149.411834717,52.9216872972
150.568775939,51.6947631836
151.23727763,50.390045166
153.265945435,49.7778711963
155.934188843,49.8835742956
158.305969238,49.5737389294
160.677734375,49.1867334409
162.675320529,48.4620666504
163.938919067,47.4491661856
163.550473712,45.841796875
161.863616943,45.0017850512

将其保存为“contour_ellipse.txt”并运行此脚本将得到下图：

import numpy
import pandas as pd
from matplotlib.patches import Ellipse

def fitEllipse(cont,method):

    x=cont[:,0]
    y=cont[:,1]

    x=x[:,None]
    y=y[:,None]

    D=numpy.hstack([x*x,x*y,y*y,x,y,numpy.ones(x.shape)])
    S=numpy.dot(D.T,D)
    C=numpy.zeros([6,6])
    C[0,2]=C[2,0]=2
    C[1,1]=-1
    E,V=numpy.linalg.eig(numpy.dot(numpy.linalg.inv(S),C))

    if method==1:
        n=numpy.argmax(numpy.abs(E))
    else:
        n=numpy.argmax(E)
    a=V[:,n]

    #-------------------Fit ellipse-------------------
    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
    cx=(c*d-b*f)/num
    cy=(a*f-b*d)/num

    angle=0.5*numpy.arctan(2*b/(a-c))*180/numpy.pi
    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)*numpy.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
    down2=(b*b-a*c)*( (a-c)*numpy.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
    a=numpy.sqrt(abs(up/down1))
    b=numpy.sqrt(abs(up/down2))

    #---------------------Get path---------------------
    ell=Ellipse((cx,cy),a*2.,b*2.,angle)
    ell_coord=ell.get_verts()

    params=[cx,cy,a,b,angle]

    return params,ell_coord

def plotConts(contour_list):
    '''Plot a list of contours'''
    import matplotlib.pyplot as plt
    fig=plt.figure()
    ax2=fig.add_subplot(111)
    for ii,cii in enumerate(contour_list):
        x=cii[:,0]
        y=cii[:,1]
        ax2.plot(x,y,'-')
    plt.show(block=False)

#-------------------Read in data-------------------
df=pd.read_csv('contour_ellipse.txt')
data=numpy.array(df)

params1,ell1=fitEllipse(data,1)
params2,ell2=fitEllipse(data,2)

plotConts([data,ell1,ell2])

小椭圆为原码，绿色为固定码。

这个错误不会每次都出现，因为很多时候最大值也是绝对值的最大值。

一些令人困惑的事情：

如果你从这里看 Fitzgibbon 的论文：http://cseweb.ucsd.edu/~mdailey/Face-Coord/ellipse-specific-fitting.pdf，他们说

由于 C 的特征值为 -2, -1, 2, 0, 0, 0，从引理 1 我们 有那个 (8) 正好有一个正特征值\lambda_i < 0

我认为这是一个错字，< 应该是 >。

另一篇讨论这个问题的论文 (https://github.com/bdhammel/least-squares-ellipse-fitting/blob/master/media/WSCG98.pdf)，他们说

我们正在寻找对应于最小值的特征向量a_k 正特征值\labmda_k

这令人困惑，因为应该只有一个正面。

最后，如果要拟合的数据数量少于 6 个，参数数量也会给你带来问题。

【参考方案2】：

我建议使用 scikit-image。

它有一个椭圆拟合函数EllipseModel，它实现了Halir, R.; Flusser, J. “椭圆的数值稳定直接最小二乘拟合”。在过程中。第六届中欧计算机图形与可视化国际会议。 WSCG（第 98 卷，第 125-132 页）.

使用您在上面定义的points 的示例：

import numpy as np
from skimage.measure import EllipseModel
from matplotlib.patches import Ellipse
import matplotlib.pyplot as plt

points = [(0,3),(1,2),(1,7),(2,2),(2,4),(2,5),(2,6),(2,14),(3,4),(4,4),(5,5),(5,14),(6,4),(7,3),(7,7),(8,10),(9,1),(9,8),(9,9),(10,1),(10,2),(10,12),(11,0),(11, 7),(12,7),(12,11),(12,12),(13,6),(13,8),(13,12),(14,4),(14,5),(14,10),(14,13)]

a_points = np.array(points)
x = a_points[:, 0]
y = a_points[:, 1]

ell = EllipseModel()
ell.estimate(a_points)

xc, yc, a, b, theta = ell.params

print("center = ",  (xc, yc))
print("angle of rotation = ",  theta)
print("axes = ", (a,b))

fig, axs = plt.subplots(2, 1, sharex=True, sharey=True)
axs[0].scatter(x,y)

axs[1].scatter(x, y)
axs[1].scatter(xc, yc, color='red', s=100)
axs[1].set_xlim(x.min(), x.max())
axs[1].set_ylim(y.min(), y.max())

ell_patch = Ellipse((xc, yc), 2*a, 2*b, theta*180/np.pi, edgecolor='red', facecolor='none')

axs[1].add_patch(ell_patch)
plt.show()

返回：

center =  (7.290242506300351, 7.317565035114109)
angle of rotation =  2.1516086051307814
axes =  (5.956065316845365, 6.195071281072721)

【讨论】：

嘿，谢谢你的回答，如果 mu 数据不是点，我如何在我的数据上使用它来绘制，如在 a_points =[...] 中定义的这里，但是两个数组 x_1（数组包含x 坐标）和 x_2（包含 y 坐标的数组）？

@Parisia，我认为您可以执行类似 a_points = np.hstack([x_1, x_2]) 的操作并从那里继续