projection orthogonale avec numpy

J'ai une liste de 3D-points pour lesquels j'ai calculer 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)

Cela me renvoie une liste de points qui sont tous dans un avion, mais quand je les affiche, ils ne sont pas à la position qu'ils devraient être.
Je crois qu'il existe déjà un certain nombre de méthode pour résoudre ce problème, mais je ne pouvais pas trouver tout =(

J'ai trouvé mon erreur: j'ai fait une fausse hypothèse: Mon calcul pour P_new est faux. C'est correct: P_new = r0 + (((x-r0)*u)/(uu))*u + (((x-r0)*w)/(ww))*w

OriginalL'auteur Munchkin | 2013-07-24

10

Vous faites une mauvaise utilisation de np.lstsq, puisque vous le nourrissez un précalculées matrice de 3x3, au lieu de le laisser faire le travail. Je voudrais faire comme ceci:
```
import numpy as np
def calc_plane(x, y, z):
a = np.column_stack((x, y, np.ones_like(x)))
return np.linalg.lstsq(a, z)[0]
>>> x = np.random.rand(1000)
>>> y = np.random.rand(1000)
>>> z = 4*x + 5*y + 7 + np.random.rand(1000)*.1
>>> calc_plane(x, y, z)
array([ 3.99795126,  5.00233364,  7.05007326])
```
Il est effectivement plus pratique d'utiliser une formule pour votre avion, qui ne dépend pas du coefficient de z n'étant pas nulle, c'est à dire utiliser a*x + b*y + c*z = 1. Vous pouvez même calculer a, b et c faire:
```
def calc_plane_bis(x, y, z):
a = np.column_stack((x, y, z))
return np.linalg.lstsq(a, np.ones_like(x))[0]
>>> calc_plane_bis(x, y, z)
array([-0.56732299, -0.70949543,  0.14185393])
```
À projeter des points sur un plan, à l'aide de mon alternative équation, le vecteur (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 fait par le référencement de tous points à ce point sur le plan, la projection des points sur la normale vecteur, soustraire que la projection de points, puis de référencement de retour à l'origine. Vous pouvez le faire comme suit:
```
def project_points(x, y, z, a, b, c):
"""
Projects the points with coordinates x, y, z onto the plane
defined by a*x + b*y + c*z = 1
"""
vector_norm = a*a + b*b + c*c
normal_vector = np.array([a, b, c]) / np.sqrt(vector_norm)
point_in_plane = np.array([a, b, c]) / vector_norm
points = np.column_stack((x, y, z))
points_from_point_in_plane = points - point_in_plane
proj_onto_normal_vector = np.dot(points_from_point_in_plane,
normal_vector)
proj_onto_plane = (points_from_point_in_plane -
proj_onto_normal_vector[:, None]*normal_vector)
return point_in_plane + proj_onto_plane
```
Alors maintenant, vous pouvez faire quelque chose comme:
```
>>> project_points(x, y, z, *calc_plane_bis(x, y, z))
array([[  0.13138012,   0.76009389,  11.37555123],
[  0.71096929,   0.68711773,  13.32843506],
[  0.14889398,   0.74404116,  11.36534936],
..., 
[  0.85975642,   0.4827624 ,  12.90197969],
[  0.48364383,   0.2963717 ,  10.46636903],
[  0.81596472,   0.45273681,  12.57679188]])
```
merci beaucoup, c'est génial!
Je viens de réaliser que le code ci-dessus est en supposant que votre plan ne passant pas par l'origine, c'est à dire qu'il ne peut pas être utilisé pour projeter sur un plan d'équation a*x + b*y + c*z = 0. Pas sûr de la façon de contourner facilement sans trop de complications, mais il faut être conscient de cette mise en garde importante.
oh merci, hier j'ai testé uniquement avec un plan qui passe par l'origine, mais vous avez raison: il n'est pas correct pour les avions en dehors de l'origine. j'ai fait la suivante:
calculer la moyenne arithmétique de tous les x,y,z-valeurs; déplacer tous les points de cette valeur (de sorte que le centre de cette scatter-plot passe par l'origine) et puis startet votre code. après que j'ai déménagé 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. J' pense il est juste à la bonne position en x, mais pas à la bonne y et z de la position de la co-système. A-il quelque chose à faire avec cette ligne: return np.linalg.lstsq(a, np.ones_like(x))[0] ?

OriginalL'auteur Jaime
1

Vous pouvez tout simplement faire tout ce qui dans des matrices est une option.

Si vous ajoutez vos points de la ligne de vecteurs à une matrice X, et y est un vecteur, alors les paramètres du vecteur beta pour la solution des moindres carrés sont:
```
import numpy as np
beta = np.linalg.inv(X.T.dot(X)).dot(X.T.dot(y))
```
mais il y a un moyen plus facile, si nous voulons faire des projections: décomposition QR nous donne un repère orthonormé matrice de projection, comme Q.T, et Q est lui-même la matrice orthonormale de vecteurs de base. Donc, on peut d'abord formulaire QR, puis obtenir beta, puis utilisez Q.T à projet de les points.

QR:
```
Q, R = np.linalg.qr(X)
```
beta:
```
# use R to solve for beta
# R is upper triangular, so can use triangular solver:
beta = scipy.solve_triangular(R, Q.T.dot(y))
```
Alors maintenant, nous avons beta, et nous pouvons projeter les points à l'aide de Q.T très simplement:
```
X_proj = Q.T.dot(X)
```
Thats it!

Si vous voulez plus d'informations et graphiques piccies et d'autres choses, j'ai fait tout un tas de notes, tout en faisant quelque chose de similaire à: https://github.com/hughperkins/selfstudy-IBP/blob/9dedfbb93f4320ac1bfef60db089ae0dba5e79f6/test_bases.ipynb

(Edit: notez que si vous souhaitez ajouter un biais terme, de sorte que le meilleur ajustement de ne pas avoir à passer par l'origine, vous pouvez simplement ajouter une colonne supplémentaire, avec tous-1s, à X, qui agit comme le biais terme/long métrage)

OriginalL'auteur Hugh Perkins

Vous devez vous connecter pour publier un commentaire.