numpyでの正射影

Question

私はnumpy.linalg.lstsq - メソッドで平面を計算する3D点のリストを持っています。しかし、今、私はこの飛行機に各点の直交射影をしたいが、私は私のミスを見つけることができません。

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はあまり使用していません。 。代わりに、それは仕事をさせるの私はこのようにそれを行うだろう：

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

は（編集：ので、あなたはバイアス項を追加する場合に注意最も適合するものは原点を通過する必要はありません。バイアス項/機能として機能するXに、すべて1の列を追加するだけです）