10

我有一个 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)

这会返回一个点列表，这些点都在一个平面上，但是当我显示它们时，它们不在应有的位置。我相信已经有某种内置方法可以解决这个问题，但是我找不到任何 =(

4

3 回答 3

14

您对 的使用非常差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])

实际上，为您的飞机使用不依赖于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])

要将点投影到平面上，使用我的替代方程，矢量(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

因此，现在您可以执行以下操作：

>>> 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]])
于 2013-07-24T15:25:05.917 回答
5

您可以简单地在矩阵中做所有事情是一种选择。

如果将点作为行向量添加到矩阵X，并且y是向量，则beta最小二乘解的参数向量为：

import numpy as np

beta = np.linalg.inv(X.T.dot(X)).dot(X.T.dot(y))

但如果我们想要进行投影，还有一种更简单的方法：QR 分解为我们提供了一个正交投影矩阵，as Q.T，并且Q它本身就是正交基向量的矩阵。所以，我们可以先 form QR，然后 get beta，然后用Q.T来投影点。

二维码：

Q, R = np.linalg.qr(X)

测试版：

# use R to solve for beta
# R is upper triangular, so can use triangular solver:
beta = scipy.solve_triangular(R, Q.T.dot(y))

所以现在我们有了beta，我们可以使用Q.T非常简单的方法来投影点：

X_proj = Q.T.dot(X)

而已！

如果你想要更多信息和图形图片和东西，我做了一大堆笔记，同时做类似的事情，在：https ://github.com/hughperkins/selfstudy-IBP/blob/9defbb93f4320ac1bfef60db089ae0dba5e79f6/test_bases.ipynb

（编辑：请注意，如果要添加偏差项，因此最佳拟合不必通过原点，您可以简单地添加一个附加列，全为 1，X作为偏差项/特征)

于 2017-02-05T23:54:33.357 回答
0

这个网页有一个非常棒的代码库。它很好地实现了Maple 所阐述的理论，具体如下：numpy

# 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)
于 2020-07-22T19:24:51.903 回答