如何使用python中的最小二乘拟合找到圆心？

Question

我正在尝试拟合一些数据点以找到圆心。以下所有点都是圆周周围的噪声数据点：

data = [(2.2176383052987667, 4.218574252410221),
(3.3041214516913033, 5.223500807396272),
(4.280815855023374, 6.461487709813785),
(4.946375258539319, 7.606952538212697),
(5.382428804463699, 9.045717060494576),
(5.752578028217334, 10.613667377465823),
(5.547729017414035, 11.92662513852466),
(5.260208374620305, 13.57722448066025),
(4.642126672822957, 14.88238955729078),
(3.820310290976751, 16.10605425390148),
(2.8099420132544024, 17.225880123445773),
(1.5731539516426183, 18.17052077121059),
(0.31752822350872545, 18.75261434891438),
(-1.2408437559671106, 19.119355580780265),
(-2.680901948575409, 19.15018791257732),
(-4.190406775175328, 19.001321726517297),
(-5.533990404926917, 18.64857428377178),
(-6.903383826792998, 17.730112542165955),
(-8.082883753215347, 16.928080323602334),
(-9.138397388219254, 15.84088004983959),
(-9.92610373064812, 14.380575762984085),
(-10.358670204629814, 13.018017342781242),
(-10.600053524240247, 11.387283417089911),
(-10.463673966507077, 10.107554951600699),
(-10.179820255235496, 8.429558128401448),
(-9.572153386953028, 7.1976672709797676),
(-8.641475289758178, 5.8312286526738175),
(-7.665976739804268, 4.782663065707469),
(-6.493033077746997, 3.8549965442534684),
(-5.092340806635571, 3.384419909199452),
(-3.6530364510489073, 2.992272643733981),
(-2.1522365767310796, 3.020780664301393),
(-0.6855406924835704, 3.0767643753777447),
(0.7848958776292426, 3.6196842530995332),
(2.0614188482646947, 4.32795711960546),
(3.2705467984691508, 5.295836809444288),
(4.359297538484424, 6.378324784240816),
(4.981264502955681, 7.823851404553242)]

我正在尝试使用像 Scipy 这样的库，但我在使用可用函数时遇到了问题。

有例如：

#  == METHOD 2 ==
from scipy      import optimize

method_2 = "leastsq"

def calc_R(xc, yc):
    """ calculate the distance of each 2D points from the center (xc, yc) """
    return sqrt((x-xc)**2 + (y-yc)**2)

def f_2(c):
    """ calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc) """
    Ri = calc_R(*c)
    return Ri - Ri.mean()

center_estimate = x_m, y_m
center_2, ier = optimize.leastsq(f_2, center_estimate)

xc_2, yc_2 = center_2
Ri_2       = calc_R(*center_2)
R_2        = Ri_2.mean()
residu_2   = sum((Ri_2 - R_2)**2)

但这似乎是在使用单个 xy？关于如何将此函数插入​​我的数据示例的任何想法？

Answer 1

您的数据点看起来相当干净，我没有看到异常值，因此许多圆拟合算法都可以工作。

我建议您从 Coope 方法开始，该方法通过神奇地线性化问题来工作：

(X-Xc)² + (Y-Yc)² = R²被改写为

2 Xc X + 2 Yc Y + R² - Xc² - Yc² = X² + Y²，然后

A X + B Y + C = X² + Y²，通过线性最小二乘法求解。

Answer 2

作为 Bas Swinckels 帖子的后续，我想我会发布我的代码来实现拟合椭圆的 Halir 和 Flusser 方法

https://github.com/bdhammel/least-squares-ellipse-fitting

使用上面的代码，你可以通过以下方法找到中心。

from ellipses import LSqEllipse
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

lsqe = LSqEllipse()
lsqe.fit(data)
center, width, height, phi = lsqe.parameters()

plt.close('all')
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.axis('equal')
ax.plot(data[0], data[1], 'ro', label='test data', zorder=1)

ellipse = Ellipse(xy=center, width=2*width, height=2*height, angle=np.rad2deg(phi),
               edgecolor='b', fc='None', lw=2, label='Fit', zorder = 2)
ax.add_patch(ellipse)

plt.legend()
plt.show()

Answer 3

我没有任何拟合圆的经验，但我使用过拟合椭圆的更一般情况。使用嘈杂的数据以正确的方式执行此操作并非易事。对于这个问题，Halir 和 Flusser 在Numerically stable direct least squares fitting of ellipses 中描述的算法效果很好。该论文包含 Matlab 代码，应该可以直接转换为 Numpy。也许你可以使用这个算法来拟合一个椭圆，然后将两个轴的平均值作为半径左右。论文中的一些参考资料也提到了拟合圆，您可能需要查找这些。

Answer 4

我知道这是一个老问题，但在 2019 年，python 中有一个圆形拟合库，名为 circle-fit。

pip install circle-fit

您可以使用两种算法之一来求解，least_squares_circle 或 hyper_fit。

import circle_fit as cf
xc,yc,r,_ = cf.least_squares_circle((data)

然后你得到xc, yc 作为解圆中心的坐标对。

Answer 5

Ian Coope 的算法（论文here）使用变量替换将问题线性化。这是迄今为止我遇到的最普遍的稳健和最快的圆形拟合算法。我已经在scikit-guess 中用python 实现了该算法。使用函数skg.nsphere_fit 进行单行解决：

>>> r, c = skg.nsphere_fit(data)

>>> r
8.138962707494084
>>> c
array([-2.45595128, 11.04352796])

这是结果图：

t = np.linspace(0, 2 * np.pi, 1000, endpoint=True)
plt.scatter(*np.array(data).T, color='r')
plt.plot(r * np.cos(t) + c[0], r * np.sin(t) + c[1])
plt.axis('equal')
plt.show()