使用Numpy查找3D点阵列中的所有4个共面点集-Java 学习之路

-1

假设我有一个 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的解决方案 .

Edit ：正如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

2 回答

以下可能不是一个非常快速的解决方案，但它的工作原理和数学/几何意义 .
但首先 - 请注意，由于"diagonal"飞机穿过你的立方体，你的例子有4个共面点的 12 子集，而不是8个 . 这可以正式化，但应该清楚（如果不是通过评论，请告诉我） .
在我们的方式中，最简单的方法是生成大小为4的所有子集（没有用于重新排序的重复），然后检查由4个点定义的体积是否为0;即这4个点中的任何3个定义包含第4个的平面 . （此方法在许多堆栈交换问题中进行了解释，并且也显示在_47508中） .

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

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组共面点 .

使用以下简单代码获得的点的简单可视化（从最后一个片段继续，带有额外的导入）：

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')

这产生了以下可视化（再次，对于此特定示例）：

希望有所帮助 .
祝好运！

回复于 2024-04-27T06:53:01+08:00

1
有四个点的70个子集，您需要计算它们形成的四面体的体积 . 如果你的形状足够接近立方体，那么共面的集合将是具有最小体积的十二个 .

对于任意体积，您还可以比较通过将体积除以四个中最大面部的面积而获得的高度 . 这将需要
```
n.(n-1).(n-2).(n-3) / 4!
```
体积计算和面积计算的四倍 .

一个详尽的方法将是可怕的（O（n ^ 4）！） . 并且矢量化将需要在几何计算之前准备顶点的所有组合 .
回复于 2024-04-27T06:53:01+08:00

使用Numpy查找3D点阵列中的所有4个共面点集

2 回答

相关问题