Wie findet man den Mittelpunkt des Kreises mit der kleinsten Python-Passung?

Question

Ich versuche, einige Datenpunkte anzupassen, um das Zentrum eines Kreises zu finden. Alle folgenden Punkte sind laut Datenpunkte um den Umfang des Kreises:

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

Ich habe versucht, einige Bibliothek wie Scipy http://wiki.scipy.org/Cookbook/Least_Squares_Circle zu verwenden, aber ich habe Probleme mit den verfügbaren Funktionen.

Es gibt zum Beispiel:

# == 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) 

Aber das scheint eine einzige xy verwendet werden? Irgendwelche Ideen, wie Sie diese Funktion an mein Datenbeispiel anschließen können?

Answer 1

Ihre Datenpunkte scheinen ziemlich sauber zu sein und ich sehe keine Ausreißer, so dass viele Kreisanpassungsalgorithmen funktionieren.

Ich empfehle Ihnen mit dem Coope Methode zu starten, die durch magische Linearisieren das Problem funktioniert:

(X-Xc)^2+(Y-Yc)^2=R² als

2XcX+2YcY+R²-Xc²-Yc²=X²+Y² neu geschrieben wird, dann

AX+BY+C=X²+Y², gelöst durch lineare kleinsten Quadrate .

Answer 2

Ich habe keine Erfahrung passende Kreise, aber ich habe mit dem allgemeineren Fall der Anpassung Ellipsen gearbeitet. Dies richtig mit verrauschten Daten zu tun ist nicht trivial. Für dieses Problem funktioniert der in Numerically stable direct least squares fitting of ellipses beschriebene Algorithmus von Halir und Flusser ziemlich gut. Das Papier enthält Matlab-Code, der einfach in Numpy übersetzt werden sollte. Vielleicht könnten Sie diesen Algorithmus verwenden, um eine Ellipse anzupassen und dann den Durchschnitt der beiden Achsen als den Radius oder so zu nehmen. Einige der Referenzen in dem Papier erwähnen auch passende Kreise, Sie könnten diese nachschlagen wollen.

Answer 3

Im Nachgang zu Bas Swinckels Post, dachte ich, ich meinen Code schreiben würde die Halir und Flusser Verfahren zum Anbringen einer Ellipse

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

Mit dem obigen Code Implementierung Sie das Zentrum finden mit die folgende Methode.

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