正交投影與numpy

Question

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

這將返回我這些都是在一個平面上點的列表，但是當我顯示他們，他們不在他們應該的位置。 我相信已經有一種內置的方法來解決這個問題，但我找不到任何=（

Answer 1

您正在做一個很差的使用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]]) 

Answer 2

你可以簡單地做一切矩陣是一個選項。

如果你爲行向量添加您的點矩陣X，並且y是一個向量，則最小二乘解的參數向量beta是：

import numpy as np 

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

，但有一個更簡單的方法，如果我們想要做投影：QR分解給我們一個正交投影矩陣，因爲Q.T和Q本身就是正交基向量的矩陣。所以，我們可以先形成QR，然後得到beta，然後用Q.T來投射點。

QR：

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) 

完蛋了！

如果您想了解更多信息和圖形piccies之類的東西，我做了一大堆筆記，同時在做類似的事情，在：https://github.com/hughperkins/selfstudy-IBP/blob/9dedfbb93f4320ac1bfef60db089ae0dba5e79f6/test_bases.ipynb

（編輯：請注意，如果你想添加一個偏置項，所以最合適的不必通過原點，你可以簡單地添加一個額外的列，所有1，到X，這是作爲偏向術語/特徵）