# Finding maximum coplanar points

I have got a set of points and I need to find the plane that would contain the maximum number of these points. Is there a simple way to do this, without trying out all permutations?

Update: I am giving an illustration below:

The red dots are the points that exist in 3d space. There can be several planes like A, B, C, D that would contain 3 or more of these points. However, the plane A contains the maximum number of them. I need to find this plane (or to be more precise, the points on this plane).

Update2: I am giving below the code I had tried earlier. This coincides with the approach suggested by @hatinacat2000 and looks to be the best available option:

import bpy
from mathutils import Vector
from itertools import combinations

def main():
o = bpy.context.object

# Take sample test points from Bezier curves
# (you need to create curves having coplanar points
# by subdividing 2d segments for example)
bps = [bp for s in o.data.splines for bp in s.bezier_points]
loc = [bp.co.copy().freeze() for bp in bps]
l = list(combinations(loc, 3))

# Map of {normal vector: 3 point combination for that normal}
# Normal rounded to 3 decimal places
normalsCos = [[Vector((round(((c[1] - c[0]).cross((c[2] - c[0])).normalized())[i], 3) \
for i in range(0,3))).freeze(), c] for c in l]

normals = [n[0] for n in normalsCos]

# Map of {count: normal}
nmap = {normals.count(n): n for n in normals}
maxCnt = max(nmap.keys())
normal = nmap[maxCnt]
coTriplets = []
for nc in normalsCos:
if(nc[0] == normal):
coTriplets.append(nc[1])

# These are the required points
planeCos = set([c for t in coTriplets for c in t])

# And the remaining points...
remainder = set(loc) - planeCos


Update3: Please see the answer from @lemon. It has got a huge performance improvement over this code (this is more of a POC), and also considers parallel planes.

• What you have is no easy problem. It is called Least-Squares Fit and I'd suggest using some already made library to solve it, like for example the Ceres solver: ceres-solver.org/index.html. Oct 2 '19 at 12:38
• Also I don't think finding exactly co-planar is what you need, because in a set there could be only 4 co-planar points giving you some random plane, while the rest of the points would sit almost on a plane but not quite. Oct 2 '19 at 12:44
• @JaroslavJerrynoNovotny Thank you for the input! Maybe I was not clear enough in my question. I am not looking for a plane that will be nearest to all the points. I don't know the math, but I guess that's what Least-Squares fit will give (and I may be totally wrong). To illustrate my problem, let's say there are 7 points, and there exists a plane - A - that contains 5 of them. There would of course be other planes that contain fewer than 5. I just want some algorithmic way to find this plane A. Oct 2 '19 at 13:20
• Are there any additional constraints, for example the distribution of the points or distance between the points on the plane? Oct 6 '19 at 16:52
• @rjg The coordinates of the points in 3d space are available, so the distance between them can be calculated. There are no constraints, these points are just randomaly distributed in space. And there is no plane as such, just the points. I want to find out the plane that has the maximum number of them. Or in other words, from the given set of points I need to find out the maximum of them that are coplanar. (In the worst case there will just be 3 such points.) Oct 7 '19 at 5:03

There is no "simple" way to do anything in numerical analysis. However, the tasks can be laid out straightforwardly.

EDIT: This idea is a brain-fart because it doesn't distinguish parallel planes. The Internet must never forget though. Pile on, people.

How about this: generate an array of normalized normal vectors (not vectors in the dynamic array sense; in the math sense. A 1x3 array representation of a floating-point ordered triple is fine. I guess you can make an object but why bother) for every 3-point combination (not permutation) of the n points you have in your set. This list will contain nC3 normal "vectors".

If no plane that can be made contains more than 3 points in the set, then your list of normal vectors will entirely contain nC3 unique normals.

If any plane exists with 4 vertices, you will have 4C3 = 4 identical normal vectors in that list.

If any plane exists with 5 vertices, you will have 5C3 = 10 identical normals.

If any plane exists with k vertices, where 3 < k < n, you will have kC3 identical normals.

You would just need to quantify the normals that are the same, identify the largest value, and deduce what k was. EDIT: Thanks to Rich Sedman for pointing out comparison should be done with dot-product. Weeding out parallel planes is another task that would need to be done optimally in the worst case scenario (a stack of parallel planes).

Doing all this would require knowing how to do a 3D cross product, how to normalize a vector, how to find a dot-product, and how to solve this equation:

k!/(3!(k-3)!) = count, Or:

k(k-1)(k-2)(k-3)!/(3!(k-3)!) = count, Or:

k(k-1)(k-2)/6 = count,

where count is the number of duplicated normal vectors, which you will have already counted before solving this. This is just a 3rd degree polynomial.

Note also that if the set is randomly chosen, and even if the range and precision of floating point values are small, it is very unlikely that you will have more than 3 points in any plane. The list would need to be curated to have more than 3 coplanar points in order to make good use of this work.

• This wouldn't distinguish parallel planes? You'd still have to check for some shared vertices to extract coplanar triangles? Oct 7 '19 at 13:27
• A plane has 2 sides. What direction do the normals of a plane point at? Oct 7 '19 at 13:50
• @Robin Betts crap, you're right...I went to sleep thinking I was going to win the Nobel Prize, lol Oct 7 '19 at 19:54
• I like this solution, and the direction of the normal shouldn't matter - to compare normals, simply use the dot product between the two vectors (which will determine how close they are to pointing the same direction), take the absolute value (so we flip those in opposite directions), then check whether equal to '1' (which indicates in same direction). To account for floating point precision, you could tweak this to check for being close to 1 (eg, above 0.9999) - provided all the points in space aren't too distant. Then just need to check for parallel planes. Oct 10 '19 at 9:08
• I had tried this approach earlier. I did not consider the parallel planes though. Have updated the post with the code snippet I had written. I am not sure if you need to solve the equations with factorials. Just max count of the same normals should do. I was hoping to avoid the combinations part but looks like it's unavoidabe. I am marking this as the right answer. Thank you! Oct 10 '19 at 15:14

The way you could do this is by trying for almost all combinations of three points, choose the plane which contains the most points. It will run in $$O(n^4)$$ for $$n$$ points.

• Thanks for your answer. I had thought about this brute force method. But this would be rather too expensive, that's the reason in my question I said 'without trying out all permutations'. Oct 5 '19 at 13:17
• I'm not sure maybe it can be done in $O(n^3)$. If you know there is one plane, it should be fast. Oct 5 '19 at 13:31
• "Is there a simple way to do this, without trying out all permutations?" Oct 8 '19 at 11:38
• Not clear, that's the question. Oct 8 '19 at 11:51

The good answer stays for hatinacat2000.

Though, the implementation could be optimized from the code given in the question, and stay in $$O(n^3)$$, but does not need to allocate combinations (and that can spare a lot of memory if there are many points).

It considers parallel planes as we progressively check best plane containing first point, then second point, then third, etc.

The code is commented below:

import bpy
from mathutils import Vector
from itertools import combinations
import random
import time

# For Bezier objects
class BezierPoint:
def __init__(self, spline, point):
self.spline = spline
self.point = point
self.co = Vector( point.co )

# For mesh objects
class MeshPoint:
def __init__(self, point):
self.point = point
self.co = Vector( point.co )

# Get the point from object (either Bezier or mesh)
def points_from_object( name ):
obj = bpy.data.objects.get(name)
if obj:
if obj.type == 'MESH':
return [MeshPoint(p) for p in obj.data.vertices]
elif obj.type == 'CURVE':
points = []
for spline in obj.data.splines:
points.extend( [BezierPoint(spline, p) for p in spline.bezier_points] )
return points
return []

# Round vector coordinates
def round_vector( vector, decimals ):
return Vector( (round(vector.x, decimals), round(vector.y, decimals), round(vector.z, decimals)) )

# Find the plane that contains the biggest amount of points
def find_plane( points ):
best_plane = None
# Go through all points
for a in range( len(points) - 2 ):
# Initialize possible planes containing a
planes = {}
pa = points[a]
# Go through all pair except a or previous
for b, c in combinations( range( a + 1, len(points) ), 2 ):
pb = points[b]
pc = points[c]
ab = pb.co - pa.co
ac = pc.co - pa.co
# The rounded normal of a b c
normal = round_vector( ab.cross(ac).normalized(), 4 ).freeze()
plane = planes.get( normal )
if not plane:
# If no previous plane with this normal initialize it
plane = set()
planes[normal] = plane
# Add the point to the plane
# Check which is best plane
best_plane = plane if best_plane is None or len(best_plane) < len(plane) else best_plane

return best_plane

print( '---------------' )

points = points_from_object( 'BezierCurve' )

start_time = time.time()

best_plane = find_plane(points)

elapsed_time = time.time() - start_time

if best_plane:
print( 'found', len(best_plane))
else:

print( "elapsed", elapsed_time )


Some enhancement in the code that should gain about 15 to 20% of performance (but with same complexity).

Only diff parts are indicated below and differences are commented:

def round_vector( vector, decimals ):
# returns a tuple so that we avoid to allocate a new vector
return (round(vector.x, decimals), round(vector.y, decimals), round(vector.z, decimals))

def find_plane( points, precision_decimals = 4 ):
best_plane = set()
planes = {} #One allocation for the set and the set is now a raw Python set
for a in range( len(points) - 2 ):
planes.clear() #Clear instead of allocating again
pa = points[a]
for b in range( a + 1, len(points) - 1 ): #Two loops instead of combinations
pb = points[b] #Get once for all 'c'
ab = pb.co - pa.co #Calculated once for all 'c'
for c in range( b + 1, len(points) ):
pc = points[c]
ac = pc.co - pa.co
normal = round_vector( ab.cross(ac).normalized(), precision_decimals ) #No freeze needed as we have a tuple
plane = planes.get( normal )
if not plane:
plane = {pa} #Immediately allocated raw set structure
planes[normal] = plane