Comme pour les lignes droites, vous pouvez utiliser la méthode des moindres carrés pour approcher un cercle à partir des points P1, P2, ..., Pn. En guise de méthode générale, je ferai référence à d'autres articles, et ici je trouverai un cercle qui passe toujours par le point de départ P1 et le point final Pn et dont l'erreur est faible. Essayez également d'implémenter une telle fonction en python (numpy).
Comme avec une ligne droite, trouvez la différence (résiduelle) entre la valeur théorique et la valeur mesurée pour chaque point, et trouvez la valeur qui réduit la somme des carrés.
Généralement, trois paramètres de coordonnée centrale $ (x, y) $ et de rayon $ r $ sont nécessaires pour représenter un cercle, mais le centre du cercle passant par les deux points $ P_1 et P_n $ est toujours le deux vertical de ces deux points. C'est sur la ligne de partage égal. Une fois le centre déterminé, le rayon est également déterminé, de sorte que les paramètres requis cette fois peuvent être unifiés.
Il n'est pas impossible d'afficher le centre et le rayon du cercle comme paramètres à l'aide d'un vecteur, mais cela compliquera le calcul, nous prendrons donc une autre méthode. Considérons le cas de $ P_1 (-a, 0), P_n (a, 0) $ pour le moment. A ce moment, le centre est le point $ (0, y_0) $ sur l'axe y, le rayon $ r ^ 2 = y_0 ^ 2 + a ^ 2 $, et l'équation du cercle à obtenir est
x^2+y^2-2yy_0-a^2=0
Ce sera. Trouvez ce $ y_0 $.
Le résidu pour chaque point est la différence entre la distance de chaque point au centre et le rayon, mais pour des raisons de simplicité de calcul, considérez la différence après auto-quadrature à l'avance. En d'autres termes, le résidu de chaque point est calculé en remplaçant $ (x_k, y_k) $ dans la formule ci-dessus.
x_k^2+y_k^2-2y_ky_0-a^2
est. La somme des carrés de ceci, c'est-à-dire
\sum \{ x_k^2+y_k^2-2y_ky_0-a^2 \} ^2
Trouvez $ y_0 $ qui minimise. Puisque cette formule est une fonction quadratique de $ y_0 $, $ y_0 $ est calculé en définissant la différenciation partielle de $ y_0 $ à 0.
\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}
Donc $ y_0 $ est:
y_0 = \frac{\sum x_k^2y_k+\sum y_k^3-a^2\sum y_k}{2\sum y_k^2}
Maintenant, pour appliquer cela à n'importe quel groupe de points, le point $ P_1 (x_1, y_1) $ va à $ (-a, 0) $ et le point $ P_n (x_n, y_n) $ pour le groupe de points. Doit être converti pour être déplacé vers $ (a, 0) $. Dans ce but,
① Effectuez un mouvement parallèle pour que le milieu $ P_c $ de $ P_1 $ et $ P_n $ se déplace vers l'origine. ② Faites une rotation pour que $ P_1 $ et $ P_n $ se déplacent sur l'axe des x.
2 étapes sont nécessaires. Je vais omettre les calculs détaillés, mais vous pouvez voir que la conversion suivante doit être effectuée.
\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)
cependant,
Je l'ai implémenté avec python.
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
Voyons si ça marche.
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
Ça s'est bien passé.
La méthode de conversion en coordonnées simples puis de calcul est plus facile à comprendre le programme, mais elle ne convient pas au calcul manuel. Mais personne ne calcule manuellement une étape aussi compliquée, donc ça va. peut être.
Recommended Posts