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

问题描述 投票：2回答：1

我有一个2D点（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()
python numpy matplotlib data-fitting
1个回答
4
投票

我认为代码中存在错误。

关于将椭圆拟合到一组数据点，我发现了一些其他问题（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.text'并运行此脚本将给出如下所示的图：

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])

enter image description here

小椭圆是原始代码，绿色椭圆是固定代码。

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

一些令人困惑的事情：

如果你从这里看看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，参数数量也会给你带来问题。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.