如何使用曲面分割坐标的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）