python - 如何在Python中找到点在曲面上的垂直投影

Question

我在 3d 空间(x、y 和 z)中有一堆点，想在 python 中找到它们在表面上的垂直投影。我的表面是使用以下函数由四个点创建的:

PRECISION = 1e-8    # Arbitrary zero for real-world purposes
def plane_from_points(points):
    centroid = np.mean(points, axis=0)
    _, eigenvalues, eigenvectors = np.linalg.svd(points - centroid)
    if eigenvalues[1] < PRECISION:
        raise ValueError("Points are aligned, can't define a plane")
    normal = eigenvectors[2]
    d = -np.dot(centroid, normal)
    plane = np.append(normal, d)
    thickness = eigenvalues[2]
    return plane, thickness

这些是我存储为列表的输入点，用于创建曲面:

surface_maker=[np.array([[44., 5., 25.],
                        [44., 25., 25.],
                        [59., 5., 75.],
                        [59., 25., 75.]])]

我想在创建的曲面中找到以下点的垂直投影:

projection_point=[np.array([[49.9,  5., 65.],
                        [48., 17., 69.],
                        [47., 25., 71.9],
                        [60., 5., 39.],
                        [55., 15., 42.1],
                        [57., 25., 40.1]])]

我尝试了以下代码来执行此操作，但它给了我水平投影，而我想找到垂直投影:

pls=[]
for i in surface_maker:
    i=i.tolist()
    pl, thickness= plane_from_points(i)
    pls.append(pl)
point_on_surf=[]
n_iter=1
for i in range (n_iter):
    a, b, c, d = pls[i][0], pls[i][1], pls[i][2], pls[i][3]
    def fun(row):
        row[0] = -(d + b * row[1] + c * row[2]) / a # calculates new x
        return row[0], row[1], row[2] # new x and old y and z
    to_be_projected=[copy.copy(projection_point[i])]
    new_points = np.asarray(list(map(fun, [x for point in to_be_projected for x in point])))
    point_on_surf.append(new_points)

在图中我想象出了我想要的东西。我只是对点使用了不同的颜色以使图形更具可读性。对于上面的三点，我显示了箭头来可视化投影。我的代码给了我红色箭头末端的点，但我想找到垂直于表面的投影点。事实上，我的代码只是为 projection_point 计算一个新的 x。图中绿色箭头显示了我想要的方向。我想对 projection_point 中存在的所有点执行此操作。 预先，我非常感谢任何帮助。

Answer 1

定义平面的一种方法是通过三个点 P、Q 和 R。(四个点不一定位于同一平面，但四个点可以。)或者，您可以通过以下方式中的一个点 P 定义平面:平面和法线向量n，您可以通过叉积确定该向量。

    n = (Q − P) × (R - P)

对范数向量进行归一化，以便获得长度为 1 的单位向量 u:

    u = n   ∕   | n |

通过单位法线u与平面P上的点的差向量的点积，可以得到点S到平面的距离d和小号:

    d = (S − P) · u

距离有符号:当S在u所在平面的一侧时为正，当S在另一侧时为负。 (它为零，当然，S 在平面上。)

通过从 S 中减去 d · u，可以得到 S 投影到平面上的点 S′:

    S′ = S − d · u = S − ((S − P) · u) · u

那么，现在让我们将其放入 Python 中。首先，点和向量类。 (它们不是一样吗？你可以这样看，但我发现区分它们很有用:两个点的差异是一个向量；一个点加一个向量是一个点；点积和叉积对于向量来说是有意义的，但是也不用于点。无论如何，我更喜欢拥有一个包含 x、y 和 z 成员的类，而不是空间向量的元组。)

无论如何，这里是:

class Point:
    def __init__(self, x, y, z):
        self.x = x
        self.y = y
        self.z = z
        
    def __repr__(self):
        return "Point(%g, %g, %g)" % (self.x, self.y, self.z)
        
    def __sub__(self, other):
        """P - Q"""        
        if isinstance(other, Vector):
            return Point(self.x - other.x,
                         self.y - other.y,
                         self.z - other.z)

        return Vector(self.x - other.x,
                      self.y - other.y,
                      self.z - other.z)
       

class Vector:
    def __init__(self, x, y, z):
        self.x = x
        self.y = y
        self.z = z
        
    def __repr__(self):
        return "Vector(%g, %g, %g)" % (self.x, self.y, self.z)
        
    def norm(self):
        """u / |u|"""        
        d = math.sqrt(self.x**2 + self.y**2 + self.z**2)
        
        return Vector(self.x / d, self.y / d, self.z / d)
        
    def __mul__(self, other):
        """dot product u · v or scaling x · u""" 
        if isinstance(other, Vector):        
            return (self.x * other.x
                  + self.y * other.y
                  + self.z * other.z)
            
        return Vector(self.x * other,
                      self.y * other,
                      self.z * other)
        

def cross(a, b):
    return Vector(a.y * b.z - a.z * b.y,
                  a.z * b.x - a.x * b.z,
                  a.x * b.y - a.y * b.x)

这些只是我们针对您的案例所需的操作。乘法有双重作用:两个向量相乘产生标量点积；将向量与标量相乘会产生缩放后的向量。

你的飞机，减少到三个点:

plane = (
    Point(44.0,  5.0, 25.0),
    Point(44.0, 25.0, 25.0),
    Point(59.0,  5.0, 75.0)
)

您想要表达的观点:

points = [
    Point(49.9,  5.0, 65.0),
    Point(48.0, 17.0, 69.0),
    Point(47.0, 25.0, 71.9),
    Point(60.0,  5.0, 39.0),
    Point(55.0, 15.0, 42.1),
    Point(57.0, 25.0, 40.1)
]

以及投影代码:

x = plane[1] - plane[0]
y = plane[2] - plane[0]
u = cross(x, y).norm()
    
for p in points:
    d = (p - plane[0]) * u
    pp = p - u * d
    
    print(pp)

这将产生投影到平面上的点:

Point(55.4963, 5, 63.3211)
Point(56.4404, 17, 66.4679)
Point(57.156, 25, 68.8532)
Point(49.1743, 5, 42.2477)
Point(49.6147, 15, 43.7156)
Point(49.2294, 25, 42.4312)

该平面无限大，因此投影点不一定位于表面的四个点之间。