如何使用曲面划分坐标的numpy数组

Question

我有两个 numpy 数组，它们都有一些点的坐标。一个数组具有切割平面的坐标 (cutting_surf)，另一个数组具有一些点的坐标 (points)：

cutting_surf=np.array([[3., 1., 3.], [2., 1., 1.],[3., 2., 3.],\
                       [2., 2., 1.], [3., 3., 3.], [2., 3., 1.]])
points=np.array([[1.2, 3., 3.], [4., 1., 2.], [2.2, 2., 1.], [1.2, 1.5, 3.], [2.5, 1.5, 3.],\
                 [2.9, 1., 3.], [2.9, 2., 2.9], [4., 3., 1.9], [3.2, 1.2, 3.], [2.2, 3., 3.5],\
                 [2.5, 3., 3.], [2.2, 1.5, 3.5]])

然后，我想从points 中删除一些多余的坐标。这些坐标太接近来自cutting_surf 的坐标。我使用了以下代码：

to_be_removed=points[np.where(np.min(distance.cdist(cutting_surf, points),axis=0)<0.5)[0],:]
cleaned_result= npi.difference(points, to_be_removed) # it removes that close points

在那之后，我的最终计划是将我的points 数组分成两个。事实上，我真的想用cutting_surf 切割points 数组。我的意思是我想让我的cleaned_result 成为：

[np.array([[2.5, 3., 3.],
          [2.5, 1.5, 3.],
          [1.2, 3., 3.],
          [1.2, 1.5, 3.],
          [2.2, 3., 3.5],
          [2.2, 1.5, 3.5]])
 np.array([[4. , 3. , 1.9],
           [4. , 1. , 2. ]])]

我的图显示了我拥有的坐标分布。我在这里只展示了一些简单的坐标，并根据它们的x 值对它们进行排序，我可以将数组分成两个，但实际上它要复杂得多。我认为唯一的方法是使用cutting_surf 坐标创建表面，然后分离两个簇。我尝试使用cutting_surf 的四个角创建曲面，但我不知道如何使用此曲面对我的点进行聚类：

def plane_from_points(cutting_surf):
    centroid = np.mean(cutting_surf, axis=0)
    _, eigenvalues, eigenvectors = np.linalg.svd(cutting_surf - 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

在此先感谢您的帮助和贡献。

Answer 1

如果我理解正确，您已经完成了距离过滤，所以我将重点介绍分离“干净”点的方面。

请注意：在您的示例中，六个cutting_surf 点位于两条平行线上，因此它们跨越一个平面 - 如果您的表面不是一个简单的平面，这种简单的方法将不起作用。

从几何上讲：

1. 我们确定与cutting_surf 平面正交的向量orth。
2. 对于points 中的每个点p，我们确定orth 和p 之间相对于平面上某个点的点积。点积的符号是平面的“边”：所有带正号的点都在一侧，所有带负号的点都在另一侧。 （而那些点积为零的人在平面上。）

对于第一部分，我们用两个范数向量v1 和v2 来描述cutting_surf 平面，而正交向量o 只是它们的叉积。这基本上是对another Stackoverflow answer的改编：

import numpy as np
# surface points slightly reordered
cutting_surf = np.array([[3., 1., 3.], [3., 2., 3.], [3., 3., 3.],
                         [2., 1., 1.], [2., 2., 1.], [2., 3., 1.]])

# take three points, not all on the same line
# -> two vectors for the plane
corner = cutting_surf[0]
v1 = np.subtract(cutting_surf[1], corner)
v1 = v1 / np.linalg.norm(v1)

v2 = np.subtract(cutting_surf[3], corner)
v2 = v2 / np.linalg.norm(v2)

orth = np.cross(v1, v2)

对于第二部分，获取每个“干净”点，确定其与平面上某个点的差异，并按与orth 的点积符号分组：

cleaned_points = np.array([[2.5, 3., 3.], [2.5, 1.5, 3.], [1.2, 3., 3.],
      [1.2, 1.5, 3.], [2.2, 3., 3.5], [2.2, 1.5, 3.5], [4. , 3. , 1.9],
      [4. , 1. , 2. ]])

one, other = [], []
for p in cleaned_points:
    # point relative to our corner:
    p_moved = np.subtract(p, corner)
    d = np.dot(orth, p_moved)
    # in case of doubt, check for d == 0 here
    (one if d < 0 else other).append(p)
print(one)
print(other)

双方one 和other 的输出似乎是所需的分组：

[array([4. , 3. , 1.9]), array([4., 1., 2.])]
[array([2.5, 3. , 3. ]), array([2.5, 1.5, 3. ]), array([1.2, 3. , 3. ]), array([1.2, 1.5, 3. ]), array([2.2, 3. , 3.5]), array([2.2, 1.5, 3.5])]

这是一个渲染不同元素的 3d 图的脚本：

from matplotlib import pyplot as plt
import numpy as np

# surface points slightly reordered
cutting_surf = np.array([[3., 1., 3.], [3., 2., 3.], [3., 3., 3.],
                        [2., 1., 1.], [2., 2., 1.], [2., 3., 1.]])

# take three points, not all on the same line
# -> two vectors for the plane
corner = cutting_surf[0]
v1 = np.subtract(cutting_surf[1], corner)
v1 = v1 / np.linalg.norm(v1)

v2 = np.subtract(cutting_surf[3], corner)
v2 = v2 / np.linalg.norm(v2)

orth = np.cross(v1, v2)

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlim((0.5, 3.5))
ax.set_ylim((0.5, 3.5))
ax.set_zlim((0.5, 3.5))

# plot the surface points
xs, ys, zs = ([p[i] for p in cutting_surf] for i in range(3))
ax.scatter3D(xs, ys, zs, color="grey", marker="+")

# our two plane vectors:
ax.plot([xs[0], xs[1]], [ys[0], ys[1]], [zs[0], zs[1]], color="grey", alpha=0.6)
ax.plot([xs[0], xs[3]], [ys[0], ys[3]], [zs[0], zs[3]], color="grey", alpha=0.6)
# the orthogonal vector, based on the corner point, in green
ax.plot([corner[0], corner[0] + orth[0]],
        [corner[1], corner[1] + orth[1]],
        [corner[2], corner[2] + orth[2]], color="green", alpha=0.6)

cleaned_points = np.array([[2.5, 3., 3.], [2.5, 1.5, 3.], [1.2, 3., 3.],
          [1.2, 1.5, 3.], [2.2, 3., 3.5], [2.2, 1.5, 3.5], [4. , 3. , 1.9],
          [4. , 1. , 2. ]])

for p in cleaned_points:
    # point relative to our corner:
    p_moved = np.subtract(p, corner)
    d = np.dot(orth, p_moved)
    ax.plot(*p, color="cyan" if d < 0 else "red", marker="D")
plt.show()

所以在这里您可以看到平面（灰色线）、正交向量（绿色）和分离的“干净”点（青色和红色）：

作为一个可能的好处，如果您对正交向量进行归一化，点积会为您提供点到平面的距离，因此您可以将清理步骤合并到这个步骤中（例如，参见https://mathinsight.org/distance_point_plane）。