Test for coplanar points, and comparing best planar fit of many sets of points?

Question

Perhaps best asked in three parts:

1. Given five not necessarily coplanar points (in 3 dimensions), what is a good measure of how close to coplanar they are?

2. Given another set of five not necessarily coplanar points, how can we assess which of these two sets of five points is “more coplanar”?

3. Given n sets of five not necessarily coplanar points, how can we order these sets of points from “most coplanar” to “least coplanar”?

Suggestions?

I’m working with sets of five points at a time, but will eventually need to consider more points in these sets.

Is this a well-formulated question? An algorithm would be helpful, especially if coded in Python.

Answer 1

If you calculate the singular value decomposition (SVD) of the point matrix, you get the semi-axes of an ellipsoid that approximates the shape of the point distribution.

If the point array has a point per row:

points = np.array([(0,0,0), (1,0,1), (0,1,0), (2,0,2), (1,1,1)])

centroid = np.mean(points, axis=0)
_, magnitudes, axes = np.linalg.svd(points - centroid, full_matrices=False)

axes is an array containing the three unitary axes of the ellipsoid, normal to each other. And magnitudes is a vector of non-negative numbers that represent the magninude of each axis. Both are conveniently sorted from high to low magnitude.

The axis of minimum magnitude is normal to the best fitting plane. Along with the centroid you get the plane equation:

>>> normal = axes[2]
>>> normal
array([-7.07106781e-01,  1.11022302e-16,  7.07106781e-01])
>>> centroid
array([0.8, 0.4, 0.8])

The magnitude of the last axis is a measure of coplanarity, where zero represents a perfect plane fit. In this example, the fit is perfect:

>>> dispersion = magnitudes[2]
>>> dispersion
1.5265928729081877e-16

You can use the dispersion value to rank point sets by coplanarity.

Answer 2

I would use a simpler approach:

1. Given your points, compute the best-fit plane.
2. Compute the point-plane distances.
3. Compute a metric based on the distances above to measure the quality of the fitting.

So, here is the code in Python (for the best-fit plane computation I use the scikit-spatial library):

Import libraries:

import numpy as np
import matplotlib.pyplot as plt

from skspatial.objects import Plane, Points

Load the points (for the sake of simplicity I have just hardcoded them here, but you will get them from an external source):

points = Points([[2, 0, -1], [1, 0, 0], [0, 4, -2], [6, -2, 5], [-2, -7, 5]])

Compute the best-fit plane:

best_fit_plane = Plane.best_fit(points)

Compute the point-plane distances:

distances = [best_fit_plane.distance_point(point) for point in points]

Compute a metric to measure the quality of the fitting (e.g, a RMSE):

error = np.sqrt(np.dot(distances, distances) / len(distances))

Now, going back to three questions:

1. Use the error metric.
2. Use the error as a comparison metric.
3. Use the distances.

The Points class has an are_coplanarmethod. So, before doing the computation above, check if the points are coplanar in the first place:

points.are_coplanar()

Answer 3

I would start with 3D OBB and try to use its dimensions as metrics m. Let a,b,c be the sides of the OBB then for example I would try this:

V = a*b*c             // volume
d = min(a,b,c)        // thickness
S = a*b*c/min(a,b,c)  // area
  = V/d
m = d/S               // cop-lanarity metrics
  = d/(V/d)
  = d^2/V

this will lead to m in range <0 , +inf> where 0 means co-planar however the result will be non linear and maybe you should normalize the result by dividing it with V so you can compare between different PCLs

If you want something linear you can try angle of side and diagonal chose the side with smallest and not biggest side lengths:

m = atan( min(a,b,c) / ( a*b*c / (min(a,b,c)*max(a,b,c)) ) )
    atan( min(a,b,c)^2 * max(a,b,c) / (a*b*c) )

this will lead to angle in range <0deg , 45deg> where 0deg means co-planar. If you want to have something more precise I would add also the other side angle and combine them somehow for example like this:

m0 = atan( min(a,b,c)^2 * max(a,b,c) / (a*b*c) ) 
m1 = atan( min(a,b,c) / max(a,b,c) )  
m =  0.5*(m0+m1)

If you sort the sides so a<=b<=c then you can rewrite to:

m0 = atan(a/b) 
m1 = atan(a/c)  
m =  0.5*(m0+m1)