最小二乘法拟合正椭圆

oZincO

已于 2023-01-11 15:29:23 修改

阅读量1.2k

点赞数

文章标签：最小二乘法 python 算法

于 2022-12-20 11:37:50 首次发布

本文链接：https://blog.csdn.net/oZincO/article/details/128381449

版权

最小二乘拟合椭圆的文章和代码很多，opencv也有现成的函数fitEllipse可以调用，但是都是拟合形式如下的椭圆形：

$Ax^{2}+Bxy+Cy^{2}+Dx+Ey+F=0$

由于有xy项，这样拟合出椭圆的长轴与x轴有夹角，由于项目需要，要做一个没有夹角的拟合，即形式如下的椭圆：

$Ax^{2}+Cy^{2}+Dx+Ey+F=0$

我这里把这种椭圆叫做“正椭圆”，查了很多资料并没有找到现成的函数，于是自己改造一下别人的函数自己实现一个，具体看了两篇文章，第一篇比较简单，但有点问题，拟合出来的只能保证是相同形式的圆锥曲线，不能保证是椭圆，第二篇文章则可以保证拟合出来的是一个椭圆

先看第一个简单但有缺陷的方法，只能保证拟合的是圆锥曲线不能保证是椭圆：

具体参考了这篇博文：最小二乘法椭圆拟合_77喵~的博客-CSDN博客_最小二乘法拟合椭圆

这是一个标准的最小二乘公式，把博文中的A取0即可，得到方程：

$\begin{bmatrix} \sum y_{i}^{4} & \sum x_{i}y_{i}^{2} & \sum y_{i}^{3} & \sum y_{i}^{2} \\ \sum x_{i}y_{i}^{2} & \sum x_{i}^{2} & \sum x_{i}y_{i} & \sum x_{i} \\ \sum y_{i}^{3} & \sum x_{i}y_{i} & \sum y_{i}^{2} & \sum y_{i} \\ \sum y_{i}^{2} & \sum x_{i} & \sum y_{i} & N \end{bmatrix} \begin{bmatrix} B\\ C\\ D\\ E\\ \end{bmatrix} = -\begin{bmatrix} \sum x_{i}^{2}y_{i}^{2} \\ \sum x_{i}^{3} \\ \sum x_{i}^{2}y_{i} \\ \sum x_{i}^{2} \end{bmatrix}$

代码：

def fit_regular_ellipse(x,y):
    '''
    l-s fit: x^2 + B*y^2 + C*x + D*y + E = 0
    该方法的缺点是只能拟合这个形式的圆锥曲线 不能保证一定是椭圆
    '''
    yi4Sum = np.sum(np.power(y,4))
    yi3Sum = np.sum(np.power(y,3))
    yi2Sum = np.sum(np.power(y,2))
    yiSum = np.sum(y)
    xi3Sum = np.sum(np.power(x,3))
    xi2Sum = np.sum(np.power(x,2))
    xiSum = np.sum(x)
    xiyiSum = np.sum(x*y)
    xiyi2Sum = np.sum(x*np.power(y,2))
    xi2yi2Sum = np.sum(np.power(x,2)*np.power(y,2))
    xi2yiSum = np.sum(np.power(x,2)*y)
    n = len(y)
    M = np.array([[yi4Sum, xiyi2Sum, yi3Sum, yi2Sum],
                  [xiyi2Sum,xi2Sum,xiyiSum,xiSum],
                  [yi3Sum,xiyiSum,yi2Sum,yiSum],
                  [yi2Sum,xiSum,yiSum,n]])
    S = - np.array([[xi2yi2Sum],[xi3Sum],[xi2yiSum],[xi2Sum]])
    R = np.linalg.solve(M,S)
    B, C, D, E, = R[0][0], R[1][0], R[2][0], R[3][0]
    inner1 = -B*C**2 - D**2 + 4*B*E
    inner2 = abs(1-B)
    x0 = - C / (2*B) 
    y0 = - D / (2*B)
    a = np.sqrt(-inner1/(2*B*(B + 1 - inner2)))
    b = np.sqrt(-inner1/(2*B*(B + 1 + inner2)))
    errArr = np.power(x,2) + B*np.power(y,2) + C*x + D*y + E 
    return x0, y0, a, b, errArr

再看第二个比较完备的方法，可以保证拟合的结果是一个椭圆：

具体参考了这篇博文：Fitting Ellipse拟合椭圆的若干方法分析 - 知乎

而这篇博文实际上主要参考了这篇文献：https://autotrace.sourceforge.net/WSCG98.pdf

github上有很多对这篇文章的实现：

GitHub - bdhammel/least-squares-ellipse-fitting: Fitting an Ellipse using a Least Squares method, in Python

GitHub - ndvanforeest/fit_ellipse: Fitting an ellipse through a set of points

现在我基于这些参考资料实现一个拟合正椭圆的函数，把文章中的B值取0，于是矩阵D变成：

$D=\begin{pmatrix} x_{1}^{2} & y_{1}^{2} & x_{1} & y_{1} & 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ x_{i}^{2} & y_{i}^{2} & x_{i} & y_{i} & 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ x_{N}^{2} & y_{N}^{2} & x_{N} & y_{N} & 1 \end{pmatrix}$

矩阵C变成：

$C=\begin{pmatrix} 0 & 2 & 0 & 0 & 0\\ 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 \end{pmatrix}$

其他部分都相同

代码：

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

def fit_regular_ellipse(x,y):
    '''
    拟合倾角为0的椭圆
    ax^2 + cy^2 + dx + fy + g = 0
    x和y是一维numpy.array
    '''
    x, y = x[:, np.newaxis], y[:, np.newaxis]
    D = np.hstack((x * x, y * y, x, y, np.ones_like(x)))
    S, C = np.dot(D.T, D), np.zeros([5, 5])
    C[0, 1], C[1, 0] = 2, 2
    U, s, V = svd(np.dot(inv(S), C))
    params = U[:, 0]
    b = 0
    c, d, f, g, a = params[1], params[2] / 2, params[3] / 2, params[4], params[0]
    num = b * b - a * c
    x0 = (c * d - b * f) / num
    y0 = (a * f - b * d) / num
    
    up = 2 * (a * f * f + c * d * d + g * b * b - 2 * b * d * f - a * c * g)
    sqt = np.sqrt(1 + 4 * b * b / ((a - c) * (a - c)))
    down1 = (b * b - a * c) * (
        (c - a) * sqt - (c + a)
    )
    down2 = (b * b - a * c) * (
        (a - c) * sqt - (c + a)
    )
    x_half_axis = np.sqrt(up / down1)
    y_half_axis = np.sqrt(up / down2)
    
    errArr = a*x*x + c*y*y + d*x + f*y + g

    return x0, y0, x_half_axis, y_half_axis, errArr