ортогональная проекция с numpy

Question

У меня есть список трехмерных точек, для которых я вычисляю плоскость методом 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) 

Thats это!

Если вы хотите получить больше информации и графическую piccies и прочее, я сделал целую кучу нот, в то время делать что-то подобное, по адресу: https://github.com/hughperkins/selfstudy-IBP/blob/9dedfbb93f4320ac1bfef60db089ae0dba5e79f6/test_bases.ipynb

(Edit: обратите внимание, что если вы хотите добавить термин смещения, так наилучшее соответствие не должно проходить через начало координат, вы можете просто добавить дополнительный столбец со всеми-1, X, который действует как термин/функция смещения)