J'ai une liste de points 3D pour lesquels je calcule un avion par numpy.linalg.lstsq - méthode. Mais maintenant, je veux faire une projection orthogonale pour chaque point dans cet avion, mais je ne trouve pas mon erreur:
from numpy.linalg import lstsq
def VecProduct(vek1, vek2):
return (vek1[0]*vek2[0] + vek1[1]*vek2[1] + vek1[2]*vek2[2])
def CalcPlane(x, y, z):
# x, y and z are given in lists
n = len(x)
sum_x = sum_y = sum_z = sum_xx = sum_yy = sum_xy = sum_xz = sum_yz = 0
for i in range(n):
sum_x += x[i]
sum_y += y[i]
sum_z += z[i]
sum_xx += x[i]*x[i]
sum_yy += y[i]*y[i]
sum_xy += x[i]*y[i]
sum_xz += x[i]*z[i]
sum_yz += y[i]*z[i]
M = ([sum_xx, sum_xy, sum_x], [sum_xy, sum_yy, sum_y], [sum_x, sum_y, n])
b = (sum_xz, sum_yz, sum_z)
a,b,c = lstsq(M, b)[0]
'''
z = a*x + b*y + c
a*x = z - b*y - c
x = -(b/a)*y + (1/a)*z - c/a
'''
r0 = [-c/a,
0,
0]
u = [-b/a,
1,
0]
v = [1/a,
0,
1]
xn = []
yn = []
zn = []
# orthogonalize u and v with Gram-Schmidt to get u and w
uu = VecProduct(u, u)
vu = VecProduct(v, u)
fak0 = vu/uu
erg0 = [val*fak0 for val in u]
w = [v[0]-erg0[0],
v[1]-erg0[1],
v[2]-erg0[2]]
ww = VecProduct(w, w)
# P_new = ((x*u)/(u*u))*u + ((x*w)/(w*w))*w
for i in range(len(x)):
xu = VecProduct([x[i], y[i], z[i]], u)
xw = VecProduct([x[i], y[i], z[i]], w)
fak1 = xu/uu
fak2 = xw/ww
erg1 = [val*fak1 for val in u]
erg2 = [val*fak2 for val in w]
erg = [erg1[0]+erg2[0], erg1[1]+erg2[1], erg1[2]+erg2[2]]
erg[0] += r0[0]
xn.append(erg[0])
yn.append(erg[1])
zn.append(erg[2])
return (xn,yn,zn)
3 Réponses :
Vous utilisez une très mauvaise utilisation de Il est en fait plus pratique d'utiliser une formule pour votre plan qui ne dépend pas du coefficient de à des points de projet Sur un plan, en utilisant mon équation alternative, le vecteur alors vous pouvez maintenant faire quelque chose comme: np.lstsq , car vous l'alimentant une matrice 3x3 précalculée, au lieu de le laisser faire le travail. Je le ferais comme ceci: z ne pas être zéro , c'est-à-dire utilisation a * x + b * y + c * z = 1 . Vous pouvez également calculer A , b et c faire: (a, b, c) est perpendiculaire au plan. Il est facile de vérifier que le point (A, B, C) / (A ** 2 + B ** 2 + C ** 2) est dans l'avion, de sorte que la projection peut être effectuée par Référencer tous les points à ce point de l'avion, faisant saillie les points sur le vecteur normal, soustrayez cette projection des points, puis les référencez-les à l'origine. Vous pouvez le faire comme suit:
@Munchkin Je viens de réaliser que le code ci-dessus suppose que votre plan ne passe pas à l'origine de l'origine, c'est-à-dire qu'il ne peut pas être utilisé pour projeter sur un plan avec équation A * x + b * y + c * z = 0 < / code>. Je ne sais pas comment travailler facilement autour de cela sans trop de complications, mais soyez conscient de cette importante mise en garde.
Oh merci hier, j'ai seulement testé avec un avion qui passe à l'origine, mais vous avez raison: ce n'est pas correct pour les avions à l'extérieur de l'origine. J'ai fait ce qui suit:
calculer la moyenne arithmétique de toutes les valeurs X, Y, Z; Déplacez tous les points par cette valeur (de sorte que le centre de cette parcelle de dispersion passe à travers l'origine), puis redémarrez votre code. Après cela, j'ai déplacé tous les points "en arrière". Maintenant, j'ai un avion qui regarde OK, mais c'est un mauvais endroit dans mon système de coordonnées. I pense Il s'agit simplement d'une position X correcte, mais pas à la position correcte Y- et Z dans le co-système. A-t-il quelque chose à voir avec cette ligne: retour np.linalg.lstsq (a, np.ones_like (x)) [0] ?
La manière standard de permettre des avions à travers l'origine est d'utiliser des coordonnées homogènes (par exemple, ajoutez une 4ème dimension qui est toujours 1).
Vous pouvez simplement tout faire dans Matrices est une option.
Si vous ajoutez vos points sous forme de vecteurs de ligne à une matrice mais il y a un moyen plus facile, si nous voulons faire des projections: la décomposition QR donne US une matrice de projection orthonormale, comme qr: bêta: donc nous avons donc nous avons thats il! Si vous voulez plus d'informations et de piccies graphiques, je fait tout un tas de notes, tout en faisant quelque chose de similaire, à: https: / /Github.com/HUGHERKINS/SHUGHERKINS/SHSTUDY-IBP/BLOB/9IDEFBB93F4320AC1BFEF60DB089AE0DBA5E79F6/TEST_BASS.IPYNB (Edit: Notez que si vous souhaitez ajouter un terme de biais, le meilleur ajustement n'a pas Pour passer à travers l'origine, vous pouvez simplement ajouter une colonne supplémentaire, avec All-1S, sur x et y est un vecteur, alors le vecteur de paramètres bêta pour la solution des moindres carrés est: qt et q est elle-même la matrice de vecteurs de base orthonormale. Donc, nous pouvons d'abord former qr , puis obtenir bêta , puis utiliser qt pour projeter les points. bêta , et nous pouvons projeter le points en utilisant qt tout simplement: X , qui agit en tant que terture / fonctionnalité de biais)
J'aime la méthode QR de résolution de la bêta. Cependant, je ne peux pas obtenir le X_Proj au travail. Q est une matrice NX3 (au moins avec le mode par défaut = 'réduit'), donc q.t.dot (x) est une matrice 3x3. Certainement pas la projection que nous voulions.
Ce Web Page a une assez bonne base de code. Il met en œuvre le Théorie exposé par érable dans Numpy Code> Très bien, comme suit: # import numpy to perform operations on vector
import numpy as np
# vector u
u = np.array([2, 5, 8])
# vector n: n is orthogonal vector to Plane P
n = np.array([1, 1, 7])
# Task: Project vector u on Plane P
# finding norm of the vector n
n_norm = np.sqrt(sum(n**2))
# Apply the formula as mentioned above
# for projecting a vector onto the orthogonal vector n
# find dot product using np.dot()
proj_of_u_on_n = (np.dot(u, n)/n_norm**2)*n
# subtract proj_of_u_on_n from u:
# this is the projection of u on Plane P
print("Projection of Vector u on Plane P is: ", u - proj_of_u_on_n)
J'ai trouvé mon erreur: j'ai fait une mauvaise hypothèse: mon calcul pour p_new a tort. Ceci est correct: p_new = r0 + ((x-r0) * u) / (u u) * u + (((x-r0) * w) / (w w) * w