python - 给定一组点 xi=(xi,yi) 在 Python 中拟合一个椭圆

Question

我正在从二维点 (x,y) 计算一系列索引。一个指标是短轴和长轴之间的比率。为了适应椭圆，我使用以下 post

当我运行这些函数时，最终结果看起来很奇怪，因为中心和轴长度与二维点不成比例

center =  [  560415.53298363+0.j  6368878.84576771+0.j]
angle of rotation =  (-0.0528033467597-5.55111512313e-17j)
axes =  [0.00000000-557.21553487j  6817.76933256  +0.j]

在此先感谢您的帮助

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

if __name__ == '__main__':

    points = [(560036.4495758876, 6362071.890493258),
     (560036.4495758876, 6362070.890493258),
     (560036.9495758876, 6362070.890493258),
     (560036.9495758876, 6362070.390493258),
     (560037.4495758876, 6362070.390493258),
     (560037.4495758876, 6362064.890493258),
     (560036.4495758876, 6362064.890493258),
     (560036.4495758876, 6362063.390493258),
     (560035.4495758876, 6362063.390493258),
     (560035.4495758876, 6362062.390493258),
     (560034.9495758876, 6362062.390493258),
     (560034.9495758876, 6362061.390493258),
     (560032.9495758876, 6362061.390493258),
     (560032.9495758876, 6362061.890493258),
     (560030.4495758876, 6362061.890493258),
     (560030.4495758876, 6362061.390493258),
     (560029.9495758876, 6362061.390493258),
     (560029.9495758876, 6362060.390493258),
     (560029.4495758876, 6362060.390493258),
     (560029.4495758876, 6362059.890493258),
     (560028.9495758876, 6362059.890493258),
     (560028.9495758876, 6362059.390493258),
     (560028.4495758876, 6362059.390493258),
     (560028.4495758876, 6362058.890493258),
     (560027.4495758876, 6362058.890493258),
     (560027.4495758876, 6362058.390493258),
     (560026.9495758876, 6362058.390493258),
     (560026.9495758876, 6362057.890493258),
     (560025.4495758876, 6362057.890493258),
     (560025.4495758876, 6362057.390493258),
     (560023.4495758876, 6362057.390493258),
     (560023.4495758876, 6362060.390493258),
     (560023.9495758876, 6362060.390493258),
     (560023.9495758876, 6362061.890493258),
     (560024.4495758876, 6362061.890493258),
     (560024.4495758876, 6362063.390493258),
     (560024.9495758876, 6362063.390493258),
     (560024.9495758876, 6362064.390493258),
     (560025.4495758876, 6362064.390493258),
     (560025.4495758876, 6362065.390493258),
     (560025.9495758876, 6362065.390493258),
     (560025.9495758876, 6362065.890493258),
     (560026.4495758876, 6362065.890493258),
     (560026.4495758876, 6362066.890493258),
     (560026.9495758876, 6362066.890493258),
     (560026.9495758876, 6362068.390493258),
     (560027.4495758876, 6362068.390493258),
     (560027.4495758876, 6362068.890493258),
     (560027.9495758876, 6362068.890493258),
     (560027.9495758876, 6362069.390493258),
     (560028.4495758876, 6362069.390493258),
     (560028.4495758876, 6362069.890493258),
     (560033.4495758876, 6362069.890493258),
     (560033.4495758876, 6362070.390493258),
     (560033.9495758876, 6362070.390493258),
     (560033.9495758876, 6362070.890493258),
     (560034.4495758876, 6362070.890493258),
     (560034.4495758876, 6362071.390493258),
     (560034.9495758876, 6362071.390493258),
     (560034.9495758876, 6362071.890493258),
     (560036.4495758876, 6362071.890493258)]


    a_points = np.array(points)
    x = a_points[:, 0]
    y = a_points[:, 1]
    from pylab import *
    plot(x,y)
    show()
    a = fitEllipse(x,y)
    center = ellipse_center(a)
    phi = ellipse_angle_of_rotation(a)
    axes = ellipse_axis_length(a)

    print "center = ",  center
    print "angle of rotation = ",  phi
    print "axes = ", axes

    from pylab import *
    plot(x,y)
    plot(center[0:1],center[1:], color = 'red')
    show()

每个顶点是一个 xi,y,i 点

二维点和拟合椭圆中心的图

使用 OpenCV 我得到以下结果:

import cv

PointArray2D32f = cv.CreateMat(1, len(points), cv.CV_32FC2)
for (i, (x, y)) in enumerate(points):
    PointArray2D32f[0, i] = (x, y)
    # Fits ellipse to current contour.
   (center, size, angle) = cv.FitEllipse2(PointArray2D32f)

(center, size, angle) 
((560030.625, 6362066.5),(10.480490684509277, 17.20206642150879),144.34889221191406)

Answer 1

fitEllipse 的计算返回虚假结果，因为 x 和 y 的值与值之间的差异相比非常大。例如，如果您尝试打印特征值 E，您会看到

array([  0.00000000e+00 +0.00000000e+00j,
         0.00000000e+00 +0.00000000e+00j,
         0.00000000e+00 +0.00000000e+00j,
        -1.36159790e-12 +8.15049878e-12j,
        -1.36159790e-12 -8.15049878e-12j,   1.18685632e-11 +0.00000000e+00j])

它们几乎都是零!显然这里存在某种数值上的不准确。

您可以通过将数据的均值移近零来解决此问题，以便值的大小更“正常”并且数字之间的差异变得更显着。

x = a_points[:, 0]
y = a_points[:, 1]
xmean = x.mean()
ymean = y.mean()
x = x-xmean
y = y-ymean

然后就可以成功找到圆心、phi和坐标轴，然后将圆心重新平移回(xmean, ymean):

center = ellipse_center(a)
center[0] += xmean
center[1] += ymean

---

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 -= xmean
    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 = [(560036.4495758876, 6362071.890493258),
     (560036.4495758876, 6362070.890493258),
     (560036.9495758876, 6362070.890493258),
     (560036.9495758876, 6362070.390493258),
     (560037.4495758876, 6362070.390493258),
     (560037.4495758876, 6362064.890493258),
     (560036.4495758876, 6362064.890493258),
     (560036.4495758876, 6362063.390493258),
     (560035.4495758876, 6362063.390493258),
     (560035.4495758876, 6362062.390493258),
     (560034.9495758876, 6362062.390493258),
     (560034.9495758876, 6362061.390493258),
     (560032.9495758876, 6362061.390493258),
     (560032.9495758876, 6362061.890493258),
     (560030.4495758876, 6362061.890493258),
     (560030.4495758876, 6362061.390493258),
     (560029.9495758876, 6362061.390493258),
     (560029.9495758876, 6362060.390493258),
     (560029.4495758876, 6362060.390493258),
     (560029.4495758876, 6362059.890493258),
     (560028.9495758876, 6362059.890493258),
     (560028.9495758876, 6362059.390493258),
     (560028.4495758876, 6362059.390493258),
     (560028.4495758876, 6362058.890493258),
     (560027.4495758876, 6362058.890493258),
     (560027.4495758876, 6362058.390493258),
     (560026.9495758876, 6362058.390493258),
     (560026.9495758876, 6362057.890493258),
     (560025.4495758876, 6362057.890493258),
     (560025.4495758876, 6362057.390493258),
     (560023.4495758876, 6362057.390493258),
     (560023.4495758876, 6362060.390493258),
     (560023.9495758876, 6362060.390493258),
     (560023.9495758876, 6362061.890493258),
     (560024.4495758876, 6362061.890493258),
     (560024.4495758876, 6362063.390493258),
     (560024.9495758876, 6362063.390493258),
     (560024.9495758876, 6362064.390493258),
     (560025.4495758876, 6362064.390493258),
     (560025.4495758876, 6362065.390493258),
     (560025.9495758876, 6362065.390493258),
     (560025.9495758876, 6362065.890493258),
     (560026.4495758876, 6362065.890493258),
     (560026.4495758876, 6362066.890493258),
     (560026.9495758876, 6362066.890493258),
     (560026.9495758876, 6362068.390493258),
     (560027.4495758876, 6362068.390493258),
     (560027.4495758876, 6362068.890493258),
     (560027.9495758876, 6362068.890493258),
     (560027.9495758876, 6362069.390493258),
     (560028.4495758876, 6362069.390493258),
     (560028.4495758876, 6362069.890493258),
     (560033.4495758876, 6362069.890493258),
     (560033.4495758876, 6362070.390493258),
     (560033.9495758876, 6362070.390493258),
     (560033.9495758876, 6362070.890493258),
     (560034.4495758876, 6362070.890493258),
     (560034.4495758876, 6362071.390493258),
     (560034.9495758876, 6362071.390493258),
     (560034.9495758876, 6362071.890493258),
     (560036.4495758876, 6362071.890493258)]

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