Finding maximum coplanar points

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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()
                plane.add(pa)
                planes[normal] = plane
            # Add the point to the plane
            plane.add(pb)
            plane.add(pc)
            # 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( 'not found' )

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
                plane.add(pb)
                plane.add(pc)
                if len(best_plane) < len(plane): #Avoids if else assignation
                    best_plane = plane

    return best_plane

Answer 4

If you don't want curved surfaces then points are on a plane when one of their (x,y,z) coords is the same e.g z (a horizontal plane). When a linear function f(x)=ax+b is involved the plane is rotated/sheared.