phép chiếu trực giao với numpy

Question

Tôi có một danh sách các điểm 3D mà tôi tính toán mặt phẳng theo phương thức numpy.linalg.lstsq. Nhưng Bây giờ tôi muốn làm một chiếu trực giao cho mỗi điểm vào máy bay này, nhưng tôi không thể tìm thấy sai lầm của tôi:

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) 

này trả về cho tôi một danh sách các điểm đó là tất cả trong một mặt phẳng, nhưng khi tôi hiển thị họ, họ không ở vị trí họ nên. Tôi tin rằng đã có một phương pháp tích hợp sẵn để giải quyết vấn đề này, nhưng tôi không thể tìm thấy bất kỳ = (

Answer 1

Bạn đang sử dụng rất kém np.lstsq, vì bạn đang cho nó một ma trận 3x3 được tính trước . thay vì để cho nó làm công việc tôi sẽ làm điều đó như thế này:

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]) 

nó thực sự là thuận tiện hơn để sử dụng một công thức cho máy bay của bạn mà không phụ thuộc vào hệ số z không phải là zero, tức là sử dụng a*x + b*y + c*z = 1. Bạn có thể tính toán tương tự a, b và c đang thực hiện:

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]) 

Để chiếu điểm lên mặt phẳng, sử dụng phương trình thay thế của tôi, vector (a, b, c) vuông góc với mặt phẳng. Thật dễ dàng để kiểm tra rằng điểm (a, b, c)/(a**2+b**2+c**2) nằm trên mặt phẳng, vì vậy có thể thực hiện phép chiếu bằng cách tham chiếu tất cả các điểm tới điểm đó trên mặt phẳng, chiếu các điểm đó vào vectơ thông thường, trừ phép chiếu đó khỏi các điểm, sau đó tham chiếu chúng trở lại nguồn gốc. Bạn có thể làm điều đó như sau:

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 

Vì vậy, bây giờ bạn có thể làm điều gì đó như:

>>> 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

Bạn chỉ có thể làm tất cả mọi thứ trong ma trận là một lựa chọn.

Nếu bạn thêm điểm của bạn như là vectơ hàng để một ma trận X, và y là một vector, sau đó các thông số vector beta cho các giải pháp bình phương nhỏ nhất là:

import numpy as np 

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

nhưng có một cách dễ dàng hơn, nếu chúng ta muốn làm dự đoán: phân rã QR cho chúng ta một ma trận chiếu trực giao, như Q.T, và Q chính nó là ma trận của vectơ cơ sở trực giao. Vì vậy, trước tiên chúng tôi có thể tạo thành QR, sau đó nhận beta, sau đó sử dụng Q.T để chiếu các điểm.

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)) 

Vì vậy, bây giờ chúng tôi có beta, và chúng tôi có thể trình chiếu các điểm sử dụng Q.T rất đơn giản:

X_proj = Q.T.dot(X) 

Thats nó!

Nếu bạn muốn biết thêm thông tin và piccies đồ họa và các công cụ, tôi đã thực hiện một bó toàn bộ các ghi chú, trong khi làm điều gì đó tương tự, tại địa chỉ: https://github.com/hughperkins/selfstudy-IBP/blob/9dedfbb93f4320ac1bfef60db089ae0dba5e79f6/test_bases.ipynb

(Edit: lưu ý rằng nếu bạn muốn thêm một thuật ngữ thiên vị, vì vậy phù hợp nhất không phải đi qua nguồn gốc, bạn có thể chỉ cần thêm một cột bổ sung, với tất cả-1, đến X, hoạt động như cụm từ/tính năng thiên vị)