8
$\begingroup$

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: enter image description here

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.

$\endgroup$
13
  • 2
    $\begingroup$ 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. $\endgroup$ Oct 2, 2019 at 12:38
  • 1
    $\begingroup$ 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. $\endgroup$ Oct 2, 2019 at 12:44
  • $\begingroup$ @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. $\endgroup$ Oct 2, 2019 at 13:20
  • 1
    $\begingroup$ Are there any additional constraints, for example the distribution of the points or distance between the points on the plane? $\endgroup$
    – Robert Gützkow
    Oct 6, 2019 at 16:52
  • $\begingroup$ @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.) $\endgroup$ Oct 7, 2019 at 5:03

4 Answers 4

3
+100
$\begingroup$

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.

$\endgroup$
8
  • 2
    $\begingroup$ This wouldn't distinguish parallel planes? You'd still have to check for some shared vertices to extract coplanar triangles? $\endgroup$
    – Robin Betts
    Oct 7, 2019 at 13:27
  • $\begingroup$ A plane has 2 sides. What direction do the normals of a plane point at? $\endgroup$ Oct 7, 2019 at 13:50
  • 3
    $\begingroup$ @Robin Betts crap, you're right...I went to sleep thinking I was going to win the Nobel Prize, lol $\endgroup$ Oct 7, 2019 at 19:54
  • 2
    $\begingroup$ 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. $\endgroup$ Oct 10, 2019 at 9:08
  • 2
    $\begingroup$ 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! $\endgroup$ Oct 10, 2019 at 15:14
2
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ 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'. $\endgroup$ Oct 5, 2019 at 13:17
  • $\begingroup$ I'm not sure maybe it can be done in $O(n^3)$. If you know there is one plane, it should be fast. $\endgroup$ Oct 5, 2019 at 13:31
  • $\begingroup$ "Is there a simple way to do this, without trying out all permutations?" $\endgroup$
    – Leander
    Oct 8, 2019 at 11:38
  • 1
    $\begingroup$ Not clear, that's the question. $\endgroup$ Oct 8, 2019 at 11:51
2
$\begingroup$

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
$\endgroup$
6
  • $\begingroup$ Should exist a better approach... will try to update later... $\endgroup$
    – lemon
    Oct 13, 2019 at 9:12
  • $\begingroup$ Thank you for the improved code, I ran it with 65 points and the difference of execution time between my code (about 22 secs) and yours (0.11 sec) was staggering! But of course there is a lot of scope for optimization in my code (I need not have created an entire map just to find the max count for example). It was more of a Proof of Concept of the approach. And both the approaches, I guess, are more or less the same. That is, taking 3 points at a time finding the normal; the normal that repeats the maximum number of times gives the plane. I look forward to the better approach. $\endgroup$ Oct 13, 2019 at 13:49
  • 1
    $\begingroup$ @BlenderDadaist, maybe I was presumptuous about "better approach". Though the code above can be optimized a bit with the same complexity. Will think a bit more about it before correcting anything.Though, I'd be interested to know the real amount of vertices you need in a real case. Is it around 100, or around 1000 or more? Note: there is a "big" difference between 2 codes: mine here considers parallel planes which is not the case in your implementation. $\endgroup$
    – lemon
    Oct 14, 2019 at 5:25
  • $\begingroup$ There would be around 100 points in practice. And I agree, there's a big difference between the two codes. $\endgroup$ Oct 14, 2019 at 6:37
  • 1
    $\begingroup$ I really appreciate your help. $\endgroup$ Oct 14, 2019 at 8:48
-1
$\begingroup$

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.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .