代码之家 › 专栏 › 技术社区 › tel

使用Numpy查找三维点阵列中4个共面点的所有集合

computational-geometry scipy numpy arrays python

tel · 技术社区 · 5 年前

n 存储在形状的Numpy数组中的三维点 (3, n) . 我想找到列表中所有4个点的集合,这样4个点是共面的。我该怎么做?

例如,给定一个 points

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

编辑

这是我用来生成的代码 要点 :

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 回复 | 直到 5 年前

ShlomiF 5 年前

以下可能不是一个很好的例子 快速的
4个共面点的子集,而不是8个,因为有“对角线”平面穿过立方体。这可以是形式化的,但应该是清楚的(如果不是通过评论让我知道)。
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组共面点。

通过以下简单代码获得的点的简单可视化(从上一个代码段继续,并有额外的导入):

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

这将产生以下可视化效果(同样,对于此特定示例):

希望有帮助。
祝你好运

Yves Daoust 5 年前

n.(n-1).(n-2).(n-3) / 4!

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

一个彻底的方法将是可怕的(O(n^4)!)。矢量化需要在正确进行几何计算之前准备所有顶点的组合。