Wie bei geraden Linien können Sie die Methode der kleinsten Quadrate verwenden, um einen Kreis aus den Punkten P1, P2, ..., Pn zu approximieren. Als allgemeine Methode werde ich auf andere Artikel verweisen, und hier werde ich einen Kreis finden, der immer durch den Startpunkt P1 und den Endpunkt Pn verläuft und dessen Fehler klein ist. Versuchen Sie auch, eine solche Funktion in Python (numpy) zu implementieren.
Ermitteln Sie wie bei einer geraden Linie die Differenz (Residuum) zwischen dem theoretischen Wert und dem gemessenen Wert für jeden Punkt und den Wert, der die Summe der Quadrate verringert.
Im Allgemeinen sind drei Parameter der Mittelkoordinate $ (x, y) $ und des Radius $ r $ erforderlich, um einen Kreis darzustellen, aber der Mittelpunkt des Kreises, der durch die beiden Punkte $ P_1 und P_n $ verläuft, sind immer die vertikalen zwei dieser beiden Punkte. Es ist auf der gleichen Trennlinie. Sobald das Zentrum bestimmt ist, wird auch der Radius bestimmt, so dass die diesmal zu erhaltenden Parameter vereinheitlicht werden können.
Es ist nicht unmöglich, den Mittelpunkt und den Radius des Kreises als Parameter mithilfe eines Vektors anzuzeigen, aber dies erschwert die Berechnung, sodass wir eine andere Methode anwenden. Betrachten wir vorerst den Fall von $ P_1 (-a, 0), P_n (a, 0) $. Zu diesem Zeitpunkt ist der Mittelpunkt der Punkt $ (0, y_0) $ auf der y-Achse, der Radius $ r ^ 2 = y_0 ^ 2 + a ^ 2 $ und die Gleichung des zu erhaltenden Kreises ist
x^2+y^2-2yy_0-a^2=0
Es wird sein. Finde dieses $ y_0 $.
Der Rest für jeden Punkt ist die Differenz zwischen dem Abstand von jedem Punkt zum Zentrum und dem Radius. Um die Berechnung zu vereinfachen, sollten Sie die Differenz nach dem Selbstquadrat im Voraus berücksichtigen. Mit anderen Worten, der Rest jedes Punktes wird berechnet, indem $ (x_k, y_k) $ in die obige Formel eingesetzt wird.
x_k^2+y_k^2-2y_ky_0-a^2
ist. Die Summe der Quadrate davon, das heißt
\sum \{ x_k^2+y_k^2-2y_ky_0-a^2 \} ^2
Finden Sie $ y_0 $, das minimiert wird. Da diese Formel eine quadratische Funktion von $ y_0 $ ist, wird $ y_0 $ berechnet, indem die partielle Differenzierung von $ y_0 $ auf 0 gesetzt wird.
\begin{align}
\sum -4y_k( x_k^2+y_k^2-2y_ky_0-a^2 ) &= 0 \\
\sum x_k^2y_k+\sum y_k^3-2y_0\sum y_k^2-a^2\sum y_k &= 0 \\
\end{align}
$ Y_0 $ ist also:
y_0 = \frac{\sum x_k^2y_k+\sum y_k^3-a^2\sum y_k}{2\sum y_k^2}
Um dies auf eine beliebige Punktgruppe anzuwenden, geht Punkt $ P_1 (x_1, y_1) $ zu $ (-a, 0) $ und Punkt $ P_n (x_n, y_n) $ für die Punktgruppe. Muss konvertiert werden, damit nach $ (a, 0) $ verschoben wird. zu diesem Zweck,
① Führen Sie eine parallele Bewegung durch, sodass sich der Mittelpunkt $ P_c $ von $ P_1 $ und $ P_n $ zum Ursprung bewegt. ② Drehen Sie so, dass sich $ P_1 $ und $ P_n $ auf der x-Achse bewegen.
Es sind 2 Schritte erforderlich. Ich werde die detaillierten Berechnungen weglassen, aber Sie können sehen, dass die folgende Konvertierung durchgeführt werden sollte.
\left( \begin{array}{c} x' \\ y' \end{array} \right) =
-\frac{1}{a}\left(
\begin{array}{cc}
x_1-x_c & -(y_1-y_c)\\
y_1-y_c & x_1-x_c
\end{array}
\right)
\left( \begin{array}{c} x-x_c \\ y-y_c \end{array} \right)
Jedoch,
Ich habe es mit Python implementiert.
def approxByArc(points):
start = points[0]
end = points[-1]
sample = points[1:-1]
midpoint = (start + end) /2
start2 = start - midpoint
a = np.linalg.norm(start2)
x = start2[0]
y = start2[1]
rotateMatrix = np.array([[x,y],[-y,x]]) * -1 / a
def conversion(point):
point = point - midpoint
p = np.array([[point[0]], [point[1]]])
p = rotateMatrix @ p
p = np.array([p[0][0],p[1][0]])
return p
sample = np.apply_along_axis(conversion, 1, sample)
xs = np.array([i[0] for i in sample])
ys = np.array([i[1] for i in sample])
p1 = np.sum(xs*xs*ys)
p2 = np.sum(ys*ys*ys)
p3 = np.sum(ys) * a * a
p4 = np.sum(ys*ys) * 2
y0 = (p1+p2-p3)/p4
center = np.array([[0],[y0]])
center = np.linalg.inv(rotateMatrix) @ center
center = np.array([center[0][0], center[1][0]])
center = center + midpoint
radius = np.linalg.norm(center - start)
return center,radius
Mal sehen, ob es funktioniert.
points = []
points.append([60,50])
n = 40
for i in range(1,n):
angle = math.pi / n * i
radius = 10 + random.uniform(-0.4,0.4)
x = 50 + radius * math.cos(angle)
y = 50 + radius * math.sin(angle)
points.append([x,y])
points.append([40,50])
r,c = approxByArc(np.array(points))
print(r)
print(c)
[50. 49.99319584]
10.00000231483068
Es ging gut.
Die Methode zum Konvertieren in einfache Koordinaten und anschließendes Berechnen ist für das Programm einfacher zu verstehen, eignet sich jedoch nicht für die manuelle Berechnung. Aber niemand berechnet einen so komplizierten Schritt manuell, also ist es okay. vielleicht.
Recommended Posts