numpy的正交投影 [英] orthogonal projection with numpy

问题描述

我有一个 3D 点列表，我通过 numpy.linalg.lstsq - 方法为其计算平面.但是现在我想对这个平面上的每个点做一个正交投影，但是我找不到我的错误:

from numpy.linalg import lstsqdef VecProduct(vek1, vek2):返回 (vek1[0]*vek2[0] + vek1[1]*vek2[1] + vek1[2]*vek2[2])def CalcPlane(x, y, z):# x, y 和 z 在列表中给出n = len(x)sum_x = sum_y = sum_z = sum_xx = sum_yy = sum_xy = sum_xz = sum_yz = 0对于范围(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 + ca*x = z - b*y - cx = -(b/a)*y + (1/a)*z - c/a'''r0 = [-c/a,0,0]u = [-b/a,1、0]v = [1/a,0,1]xn = []y = []zn = []# 用 Gram-Schmidt 正交化 u 和 v 得到 u 和 wuu = VecProduct(u, u)vu = VecProduct(v, u)fak0 = vu/uuerg0 = [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对于范围内的 i(len(x)):xu = VecProduct([x[i], y[i], z[i]], u)xw = VecProduct([x[i], y[i], z[i]], w)fak1 = xu/uufak2 = xw/wwerg1 = [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]]尔格[0] += r0[0]xn.append(erg[0])yn.append(erg[1])zn.append(erg[2])返回 (xn,yn,zn)

这会返回一个平面上所有点的列表，但是当我显示它们时，它们不在它们应该位于的位置.我相信已经有一定的内置方法可以解决这个问题，但是我找不到任何 =(

解决方案

你对 np.lstsq 的使用很糟糕，因为你给它提供了一个预先计算的 3x3 矩阵，而不是让它做这份工作.我会这样做:

将 numpy 导入为 npdef calc_plane(x, y, z):a = np.column_stack((x, y, np.ones_like(x)))返回 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)数组([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))返回 np.linalg.lstsq(a, np.ones_like(x))[0]>>>calc_plane_bis(x, y, z)数组([-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):"""将坐标为 x、y、z 的点投影到平面上由 a*x + b*y + c*z = 1 定义"""vector_norm = a*a + b*b + c*cnormal_vector = np.array([a, b, c])/np.sqrt(vector_norm)point_in_plane = np.array([a, b, c])/vector_norm点数 = np.column_stack((x, y, z))points_from_point_in_plane = 点 - point_in_planeproj_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)返回 point_in_plane + proj_onto_plane

所以现在您可以执行以下操作:

<预><代码>>>>project_points(x, y, z, *calc_plane_bis(x, y, z))数组([[ 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]])

I have a list of 3D-points for which I calculate a plane by numpy.linalg.lstsq - method. But Now I want to do a orthogonal projection for each point into this plane, but I can't find my mistake:

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)

This returns me a list of points which are all in a plane, but when I display them, they are not at the positions they should be. I believe there is already a certain built-in method to solve this problem, but I couldn't find any =(

解决方案

You are doing a very poor use of np.lstsq, since you are feeding it a precomputed 3x3 matrix, instead of letting it do the job. I would do it like this:

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

It is actually more convenient to use a formula for your plane that doesn't depend on the coefficient of z not being zero, i.e. use a*x + b*y + c*z = 1. You can similarly compute a, b and c doing:

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

To project points onto a plane, using my alternative equation, the vector (a, b, c) is perpendicular to the plane. It is easy to check that the point (a, b, c) / (a**2+b**2+c**2) is on the plane, so projection can be done by referencing all points to that point on the plane, projecting the points onto the normal vector, subtract that projection from the points, then referencing them back to the origin. You could do that as follows:

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

So now you can do something like:

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

这篇关于numpy的正交投影的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！