[PYTHON] Approximation par la méthode des moindres carrés d'un cercle à deux points fixes

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

politique

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.

Méthode

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}

Transformation de coordonnées

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,x_c,y_cEstP_cCoordonnées,aEst|P_1P_c|est. Après avoir effectué cette conversion sur tous les points, effectuez une approximation et effectuez la conversion inverse pour obtenir le centre et le rayon.

la mise en oeuvre

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é.

en conclusion

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

Approximation par la méthode des moindres carrés d'un cercle à deux points fixes
Distance approximative entre deux points à la surface d'un ellipsoïde en rotation (à la surface de la terre)
Calcul de la matrice d'homographie par la méthode des moindres carrés (méthode DLT)
Hit une méthode d'une instance de classe avec l'API Web Python Bottle
Dessiner un cercle passant par deux points
Analyse de régression simple par la méthode des moindres carrés
Essayez d'envoyer les résultats agrégés des deux enregistrements par e-mail avec pykintone
Trouvez le coefficient du polypole le moins carré
À propos de la précision de la méthode de calcul du rapport de circonférence d'Archimède
J'ai essayé la méthode des moindres carrés en Python
Approximer une courbe de Bézier à travers un point spécifié en utilisant la méthode des moindres carrés en Python
Détruire l'expression intermédiaire de la méthode sweep avec Python
Prenez des captures d'écran LCD avec Python-LEGO Mindstorms
Points Python du point de vue d'un programmeur en langage C
Visualisez le vocabulaire caractéristique d'un document avec D3.js
Calculer le produit des matrices avec une expression de caractère?
La méthode de copie de pandas.DataFrame est une copie profonde par défaut
Créez un tableau à deux dimensions en ajoutant une ligne à la fin d'un tableau vide avec numpy
Comment tracer beaucoup de légendes en changeant la couleur du graphique en continu avec matplotlib
Sauvegardez la sortie de GAN une par une ~ Avec l'implémentation de GAN par PyTorch ~
Une implémentation Python simple de la méthode k-voisinage (k-NN)
Un diagramme de réseau a été créé avec les données du COVID-19.
Mesurer l'importance des entités avec un outil de forêt aléatoire
Obtenez l'identifiant d'un GPU avec une faible utilisation de la mémoire
Obtenez UNIXTIME au début d'aujourd'hui avec une commande
Analysez le modèle thématique pour devenir romancier avec GensimPy3
Méthode d'approximation numérique lorsque le calcul de la dérivée est gênant
L'histoire de la création d'un bot de boîte à questions avec discord.py
Script Python qui compare le contenu de deux répertoires
Prouvons le théorème d'addition d'une fonction triangulaire en remplaçant la fonction par une fonction dans SymPy (≠ substitution)
Comment intercepter ou falsifier la communication SSL de l'appareil iOS réel par un proxy
[Python] Comment utiliser l'instruction for. Une méthode d'extraction en spécifiant une plage ou des conditions.