CODEWITHSUNDEEP
pythonarraysnumpyscipycomputational-geometry
tel
tel

Reputation: 13999

Find all sets of 4 coplanar points in an array of 3D points using Numpy

Say that I have a list of n 3D points stored in a Numpy array of shape (3, n). I want to find all of the sets of 4 points in that list such that the 4 points are coplanar. How can I do that?

For example, given an array of points containing the 8 vertices (in no particular order) of a cube rotated about an arbitrary angle in 3D space:

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

how do I find the 6 sets of 4 vertices that lie at the corners of each face of the cube? Specifically, I'm looking for a fully vectorized Numpy/Scipy based solution.

Edit: as ShlomiF points out, there are actually 12 coplanar sets of a cube's vertices, including the vertices that lie on the planes along the face diagonals of the cube.

Here's the code I used to generate 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

Upvotes: 1

Views: 1132

Answers (2)

user1196549
user1196549

Reputation:

There are 70 subsets of four points and you need to compute the volumes of the tetrahedra they form. If your shape is sufficiently close to a cube, the co-planar sets will be the twelve ones with the smallest volume.

For an arbitrary volume, you can also compare the heights obtained by dividing the volume by the area of the largest face among four. This will take

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

volume computations and four times as many area computations.

An exhaustive approach will be terrible (O(n^4) !). And vectorization will require to prepare all combinations of the vertices before the geometric computation proper.

Upvotes: 0

ShlomiF
ShlomiF

Reputation: 2895

The following might not be a very fast solution, but it works and makes mathematical/geometrical sense.
But first - please note that your example has 12 subsets of 4 coplanar points, not 8, due to having "diagonal" planes going through your cube. That can be formalized but should be clear as is (let me know if not via the comments).
That out of our way, the simplest method would be to generate all subsets of size 4 (without repeats for reordering), and then checking if the volume defined by the 4 points is 0; i.e. any 3 of those 4 points define a plane containing the 4th. (This method is explained in many stack-exchange questions, and also shows up in the wolfram definition of "Coplanar").

Implementing this can be done really simply as follows: 

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"

And you get all 12 4-sets of co-planar points.

A simple visualization of the points obtained with the following simple code (continued from the last snippet, with extra imports): 

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

Which produces the following visualization (again, for this specific example):

12 co-planar points

Hope that helps.
Good luck!

Upvotes: 1

Related Questions