numpy로 직교 투영

Question

나는 numpy.linalg.lstsq - method로 평면을 계산하는 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의 사용 빈도가 매우 낮습니다. 미리 계산 된 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

(편집 : 당신은 바이어스 용어를 추가 할 경우주의 때문에 가장 적합하지 않은 점은 원점을 통과해야합니다. 바이어스 용어/기능으로 작동하는 X에 all-1을 포함한 추가 열을 추가하면됩니다.