最小二乘拟合椭圆的文章和代码很多,opencv也有现成的函数fitEllipse可以调用,但是都是拟合形式如下的椭圆形:
由于有xy项,这样拟合出椭圆的长轴与x轴有夹角,由于项目需要,要做一个没有夹角的拟合,即形式如下的椭圆:
我这里把这种椭圆叫做“正椭圆”,查了很多资料并没有找到现成的函数,于是自己改造一下别人的函数自己实现一个,具体看了两篇文章,第一篇比较简单,但有点问题,拟合出来的只能保证是相同形式的圆锥曲线,不能保证是椭圆,第二篇文章则可以保证拟合出来的是一个椭圆
先看第一个简单但有缺陷的方法,只能保证拟合的是圆锥曲线不能保证是椭圆:
具体参考了这篇博文:最小二乘法椭圆拟合_77喵~的博客-CSDN博客_最小二乘法拟合椭圆
这是一个标准的最小二乘公式,把博文中的A取0即可,得到方程:
代码:
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 - ndvanforeest/fit_ellipse: Fitting an ellipse through a set of points
现在我基于这些参考资料实现一个拟合正椭圆的函数,把文章中的B值取0,于是矩阵D变成:
矩阵C变成:
其他部分都相同
代码:
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