Найти все наборы из 4 копланарных точек в массиве трехмерных точек, используя Numpy

Question

Скажите, что у меня есть список n 3D-точек, хранящихся в массиве Numpy формы (3, n).Я хочу найти все наборы из 4 точек в этом списке так, чтобы 4 точки были копланарными.Как я могу это сделать?

Например, учитывая массив points, содержащий 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

Answer 1

Существует 70 подмножеств из четырех точек, и вам нужно вычислить объемы тетраэдров, которые они образуют.Если ваша форма достаточно близка к кубу, копланарные множества будут двенадцатью с наименьшим объемом.

Для произвольного объема вы также можете сравнить высоты, полученные путем деления объема на площадьсамого большого лица среди четырех.Это займет

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

вычислений объема и в четыре раза больше вычислений площади.

Исчерпывающий подход будет ужасен (O (n ^ 4)!).А векторизация потребует подготовки всех комбинаций вершин перед собственно геометрическим вычислением.

Answer 2

Следующее может не быть очень быстрым решением, но оно работает и имеет математический / геометрический смысл.
Но сначала - обратите внимание, что в вашем примере есть 12 подмножеств4 копланарные точки, а не 8, из-за наличия «диагональных» плоскостей, проходящих через ваш куб.Это может быть формализовано, но должно быть ясно, как есть (дайте мне знать, если не через комментарии).
Это, по нашему мнению, самый простой способ будет генерировать все подмножества размера 4 (без повторов для переупорядочения), изатем проверка, равен ли объем, определенный 4 точками, 0;т.е. любые 3 из этих 4 точек определяют плоскость, содержащую 4-ю.(Этот метод объясняется во многих вопросах обмена стека, а также показан в определении «wolfram» «Копланарность» ).

Реализовать это можно на самом деле просто следующим образом:

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

, который создает следующую визуализацию (опять же, для этого конкретногопример):

Надеюсь, это поможет.
Удачи!