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

Question

假设我有一个n 3D 点列表，存储在形状为(3, n) 的 Numpy 数组中。我想在该列表中找到所有 4 个点的集合，以使 4 个点共面。我该怎么做？

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

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 个 4 个共面点的子集，而不是 8 个，因为有“对角线”平面穿过您的立方体。这可以正式化，但应该清楚（如果不是通过 cmets，请告诉我）。
对我们来说，最简单的方法是生成所有大小为 4 的子集（不重复重新排序），然后检查 4 个点定义的体积是否为 0；即这 4 个点中的任何 3 个定义了一个包含第 4 个的平面。 （此方法在许多堆栈交换问题中都有说明，也出现在 the wolfram definition of "Coplanar" 中）。

实现这一点可以很简单地完成如下：

import numpy as np
import scipy.linalg as spl
from itertools import combinations

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

subsets_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]  # due to precision stuff, you cant just do "det(thing) == 0"

你得到所有 12 个 4 组共面点。

使用以下简单代码获得的点的简单可视化（接上一个 sn-p，带有额外的导入）：

import matplotlib.pyplot as plt
from 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):

    ax3d = plt.subplot(3, 4, i+1, projection='3d')

    # Draw cube:
    for point_pair in neighbouring_points:
        ax3d.plot(point_pair[0], point_pair[1], point_pair[2], 'k')

    # Choose coplanar set:    
    p = coplanar_points[i]

    # Draw set:
    for x, y, z in p:
        ax3d.scatter(x, y, z, s=30, c='m')
    ax3d.set_xticks([])
    ax3d.set_yticks([])
    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) / 4!

体积计算和四倍的面积计算。

详尽的方法将是可怕的（O（n ^ 4）！）。并且矢量化需要在适当的几何计算之前准备好所有顶点组合。