Python numpy的正交投影
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/17836880/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
orthogonal projection with numpy
提问by Munchkin
I have a list of 3D-points for which I calculate a plane by numpy.linalg.lstsq - method. But Now I want to do a orthogonal projection for each point into this plane, but I can't find my mistake:
我有一个 3D 点列表，我通过 numpy.linalg.lstsq - 方法为其计算平面。但是现在我想对这个平面的每个点进行正交投影，但是我找不到我的错误：
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)
This returns me a list of points which are all in a plane, but when I display them, they are not at the positions they should be. I believe there is already a certain built-in method to solve this problem, but I couldn't find any =(
这会返回一个点列表，这些点都在一个平面上，但是当我显示它们时，它们不在它们应该位于的位置。我相信已经有某种内置方法可以解决这个问题，但是我找不到任何 =(
采纳答案by Jaime
You are doing a very poor use of np.lstsq, since you are feeding it a precomputed 3x3 matrix, instead of letting it do the job. I would do it like this:
您对 的使用非常糟糕np.lstsq，因为您正在为它提供一个预先计算的 3x3 矩阵，而不是让它完成工作。我会这样做：
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])
It is actually more convenient to use a formula for your plane that doesn't depend on the coefficient of znot being zero, i.e. use a*x + b*y + c*z = 1. You can similarly compute a, band cdoing:
实际上，为您的飞机使用不依赖于z非零系数的公式更方便，即使用a*x + b*y + c*z = 1. 您可以类似地计算a，b并c执行以下操作：
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])
To project points onto a plane, using my alternative equation, the vector (a, b, c)is perpendicular to the plane. It is easy to check that the point (a, b, c) / (a**2+b**2+c**2)is on the plane, so projection can be done by referencing all points to that point on the plane, projecting the points onto the normal vector, subtract that projection from the points, then referencing them back to the origin. You could do that as follows:
要将点投影到平面上，使用我的替代方程，向量(a, b, c)垂直于平面。很容易检查该点(a, b, c) / (a**2+b**2+c**2)是否在平面上，因此可以通过参考平面上该点的所有点，将这些点投影到法向量上，从这些点中减去该投影，然后将它们参考回到平面上来完成投影起源。你可以这样做：
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
So now you can do something like:
所以现在您可以执行以下操作：
>>> 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]])
回答by Hugh Perkins
You can simply do everything in matrices is one option.
您可以简单地在矩阵中完成所有操作是一种选择。
If you add your points as row vectors to a matrix X, and yis a vector, then the parameters vector betafor the least squares solution are:
如果您将点作为行向量添加到矩阵中X，并且y是向量，则beta最小二乘解的参数向量为：
import numpy as np
beta = np.linalg.inv(X.T.dot(X)).dot(X.T.dot(y))
but there's an easier way, if we want to do projections: QR decomposition gives us an orthonormal projection matrix, as Q.T, and Qis itself the matrix of orthonormal basis vectors. So, we can first form QR, then get beta, then use Q.Tto project the points.
但是有一个更简单的方法，如果我们想要做投影：QR 分解给了我们一个正交投影矩阵，as Q.T，Q它本身就是正交基向量的矩阵。所以，我们可以先形成QR，然后得到beta，然后用Q.T来投影点。
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))
So now we have beta, and we can project the points using Q.Tvery simply:
所以现在我们有了beta，我们可以Q.T非常简单地投影点：
X_proj = Q.T.dot(X)
Thats it!
就是这样！
If you want more information and graphical piccies and stuff, I made a whole bunch of notes, whilst doing something similar, at: https://github.com/hughperkins/selfstudy-IBP/blob/9dedfbb93f4320ac1bfef60db089ae0dba5e79f6/test_bases.ipynb
如果你想要更多的信息和图形图片和东西，我做了一大堆笔记，同时做类似的事情，在：https: //github.com/hughperkins/selfstudy-IBP/blob/9dedfbb93f4320ac1bfef60db089ae0dba5e79f6/test_bases.ipynb
(Edit: note that if you want to add a bias term, so the best-fit doesnt have to pass through the origin, you can simply add an additional column, with all-1s, to X, which acts as the bias term/feature)
（编辑：请注意，如果您想添加一个偏置项，因此最佳拟合不必通过原点，您可以简单地添加一个额外的列，全 1， to X，作为偏置项/特征)
