2014-10-26 844 views
5

我想要適合一些數據點,以找到一個圓的中心。以下所有的點都圍繞着圓的圓周噪聲數據點:如何使用python中的最小二乘擬合找到圓的中心?

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.22588), 
(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的http://wiki.scipy.org/Cookbook/Least_Squares_Circle,但我在使用的可用功能的麻煩。

有例如:

# == 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?有關如何將此函數插入​​到我的數據示例的任何想法?

回答

2

您的數據點看起來相當乾淨,而且我看不到任何異常值,所以很多圓圈擬合算法都可以使用。

我建議你開始與Coope方法,它的工作原理是神奇線性化問題:

(X-Xc)^2+(Y-Yc)^2=R²被改寫爲

2XcX+2YcY+R²-Xc²-Yc²=X²+Y²,然後

AX+BY+C=X²+Y²,通過線性最小二乘法解決。

2

我沒有任何經驗擬合圓,但我已經與更一般的擬合橢圓的情況。用嘈雜的數據以正確的方式做到這一點不是微不足道的。對於這個問題,Halir和Flusser在Numerically stable direct least squares fitting of ellipses中描述的算法工作得很好。該論文包括Matlab代碼,應該直接轉換成Numpy。也許你可以用這個算法來擬合一個橢圓,然後以兩個軸的平均值爲半徑左右。文章中的一些參考文獻還提到了擬合圓,您可能想要查看它們。

+0

這實在是不平凡的。看看你的建議。 – mimoralea 2014-10-26 19:32:38

3

作爲後續行動,以巴斯Swinckels後,我想我會後我的代碼實現擬合橢圓

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

使用上述代碼的Halir和Flusser方法,你可以找到與中心以下方法。

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

enter image description here