使用 Numpy 在 3D 点数组中查找所有 4 个共面点集

Question

假设我有一个n3D 点列表存储在一个 Numpy 数组中(3, n)。我想在该列表中找到所有 4 个点的集合，使 4 个点共面。我怎样才能做到这一点？

例如，给定一个points包含 8 个顶点的数组（无特定顺序）在 3D 空间中围绕任意角度旋转的立方体：

points = np.array([[ 0.8660254 ,  0.8660254 ,  0.        ,  0.3660254 , -0.5       ,  0.3660254 ,  0.        , -0.5       ],

                   [ 0.35355339, -0.35355339,  0.70710678, -0.25881905,  0.09473435, -0.96592583,  0.        , -0.61237244],

                   [ 1.06066017,  0.35355339,  0.70710678,  1.67303261,  1.31947922,  0.96592583,  0.        ,  0.61237244]])

我如何找到位于立方体每个面的角落的 6 组 4 个顶点？具体来说，我正在寻找一个完全矢量化的基于 Numpy/Scipy 的解决方案。

编辑：正如 ShlomiF 指出的那样，立方体的顶点实际上有 12 个共面组，包括沿着立方体的面对角线位于平面上的顶点。

这是我用来生成的代码points：

import numpy as np

import scipy.linalg as spl

def rot(axis, theta):

    return spl.expm(np.cross(np.eye(len(axis)), axis/spl.norm(axis)*theta))

rot3 = rot((1,0,0), np.pi/4) @ rot((0,1,0), np.pi/3) @ rot((0,0,1), np.pi/2)

points = np.array([[1, 0, 1, 1, 1, 0, 0, 0],

                   [0, 0, 0, 1, 1, 1, 0, 1],

                   [1, 1, 0, 1, 0, 1, 0, 0]])

points = rot3 @ points

Answer 1

以下可能不是一个非常快速的解决方案，但它有效并且具有数学/几何意义。但首先 - 请注意，您的示例有12个共面点的12个子集，而不是 8 个，因为“对角线”平面穿过您的立方体。这可以正式化，但应该保持原样（如果不是通过评论让我知道）。最简单的方法是生成大小为 4 的所有子集（无需重复重新排序），然后检查由 4 个点定义的体积是否为 0；即这 4 个点中的任何 3 个点定义了一个包含第 4 个点的平面。（此方法在许多堆栈交换问题中都有解释，并且也出现在"Coplanar" 的 wolfram 定义中）。实现这一点非常简单，如下所示：import numpy as npimport scipy.linalg as splfrom itertools import combinationsdef rot(axis, theta):&nbsp; &nbsp; return spl.expm(np.cross(np.eye(len(axis)), axis/spl.norm(axis)*theta))rot3 = rot((1,0,0), np.pi/4) @ rot((0,1,0), np.pi/3) @ rot((0,0,1), np.pi/2)points = np.array([[1, 0, 1, 1, 1, 0, 0, 0],&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;[0, 0, 0, 1, 1, 1, 0, 1],&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;[1, 1, 0, 1, 0, 1, 0, 0]])points = rot3 @ pointssubsets_of_4_points = list(combinations(points.T, 4)) # 70 subsets. 8 choose 4 is 70.coplanar_points = [p for p in subsets_of_4_points if np.abs(np.linalg.det(np.vstack([np.stack(p).T, np.ones((1, 4))]))) < 0.000001]&nbsp; # due to precision stuff, you cant just do "det(thing) == 0"并且您将获得所有 12 个 4 组共面点。使用以下简单代码获得的点的简单可视化（从最后一个片段继续，有额外的导入）：import matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3D# Get pairs of points for plotting the lines of the cube:all_pairs_of_points = list(combinations(points.T, 2))# Keep only points with distance equal to 1, to avoid drawing diagonals:neighbouring_points = [list(zip(list(p1), list(p2))) for p1, p2 in all_pairs_of_points if np.abs(np.sqrt(np.sum((p1 - p2)**2)) - 1) < 0.0001]plt.figure()for i in range(12):&nbsp; &nbsp; ax3d = plt.subplot(3, 4, i+1, projection='3d')&nbsp; &nbsp; # Draw cube:&nbsp; &nbsp; for point_pair in neighbouring_points:&nbsp; &nbsp; &nbsp; &nbsp; ax3d.plot(point_pair[0], point_pair[1], point_pair[2], 'k')&nbsp; &nbsp; # Choose coplanar set:&nbsp; &nbsp;&nbsp;&nbsp; &nbsp; p = coplanar_points[i]&nbsp; &nbsp; # Draw set:&nbsp; &nbsp; for x, y, z in p:&nbsp; &nbsp; &nbsp; &nbsp; ax3d.scatter(x, y, z, s=30, c='m')&nbsp; &nbsp; ax3d.set_xticks([])&nbsp; &nbsp; ax3d.set_yticks([])&nbsp; &nbsp; ax3d.set_zticks([])plt.suptitle('Coplanar sets of 4 points of the rotated 3D cube')这产生了以下可视化（同样，对于这个特定的例子）：

Answer 2

有四个点的 70 个子集，您需要计算它们形成的四面体的体积。如果您的形状足够接近立方体，则共面组将是体积最小的十二组。对于任意体积，您还可以比较通过将体积除以四个面中最大面的面积而获得的高度。这将需要n.(n-1).(n-2).(n-3)&nbsp;/&nbsp;4!体积计算和四倍的面积计算。详尽的方法会很糟糕（O（n^4）！）。并且矢量化需要在适当的几何计算之前准备顶点的所有组合。