Profile Matching in a Point Cloud

Question

A point cloud is generated using uniform random function for (x,y,z). As shown on the following figure, a flat intersecting plane (profile) is being investigated that matches as the best (even if not the exact one) a target profile i.e., given at the left-bottom corner. So the question is:

1- How to find such a match of given target 2D point map through point cloud considering the following notes/conditions?
2- What are then the coordinates/orientations/degree of similarity etc?

Note 1: The profile of interest could be anywhere with any rotation along axes and also could be in different shape e.g., a triangle, rectangle, quadrangle etc depending on its location and orientation. At the following demonstration only a simple rectangle one is shown.

Note 2: A tolerance value could be considered as the distance of the points from the profile. To demonstrate this for the following figure suppose a tolerance of 0.01 times the smallest dimension (~1) so tol=0.01. So if we remove the rest and project all remaining points on the plane of the profile being investigated then we will be able to check its similarity with the target profile.

Note 3: A related topic could be found at Point Pattern Recognition.

Answer 1

This is always going to require a lot of computation, especially if you want to process as many as 2000 points. I'm sure there are highly-optimized solutions for this kind of pattern-matching already, but you have to figure out what it's called in order to find them.

Since you're talking about a point cloud (sparse data) instead of an image, my cross-correlation method doesn't really apply (and would be even worse computationally). Something like RANSAC probably finds a match quickly, but I don't know much about it.

My attempt at a solution:

Assumptions:

• You want to find the best match, not just a loose or "probably correct" match
• Match will have some small amount of error due to noise in measurement or computation
• Source points are coplanar
• All source points must exist in target (= any unmatched point is a mismatch for the entire profile)

So you should be able to take a lot of shortcuts by disqualifying things and decrease computation time. In short:

1. pick three points from source
2. search through target points, finding sets of 3 points with the same shape
3. when a match of 3 points is found, check all the other points in the plane they define to see if they're a close match
4. if more than one match of all points is found, choose the one with the smallest sum of 3D distances error

More detailed:

pick a point from the source for testing s1 = (x1, y1)
Find nearest point in source s2 = (x2, y2)
d12 = (x1-x2)^2 + (y1-y2)^2
Find second nearest point in source s3 = (x3, y3)
d13 = (x1-x3)^2 + (y1-y3)^2
d23 = (x2-x3)^2 + (y2-y3)^2

for all (x,y,z) test points t1 in target:
    # imagine s1 and t1 are coincident
    for all other points t2 in target:
        if distance from test point > d12:    
            break out of loop and try another t2 point
        if distance ≈ d12:
            # imagine source is now rotated so that s1 and s2 are collinear with t1 and t2
            for all other points t3 in target:
                if distance from t1 > d13 or from t2 > d23:
                    break and try another t3
                if distance from t1 ≈ d13 and from t2 ≈ d23:
                    # Now you've found matching triangles in source and target
                    # align source so that s1, s2, s3 are coplanar with t1, t2, t3
                    project all source points onto this target plane 
                    for all other points in source:
                        find nearest point in target
                        measure distance from source point to target point
                        if it's not within a threshold:
                            break and try a new t3
                        else:
                            sum errors of all matched points for this configuration (defined by t1, t2, t3)

Whichever configuration has the least squared error for all other points is the best match

Since we're working with 3 nearest neighbor test points, matching target points can be simplified by checking if they are within some radius. If searching for a radius of 1 from (0, 0), for instance, we can disqualify (2, 0) based on x1 - x2, without calculating the actual Euclidean distance, to speed it up a bit. This assumes that subtraction is faster than multiplication. There are optimized searches based on more arbitrary fixed radius, too.

function is_closer_than(x1, y1, z1, x2, y2, z2, distance):
    if abs(x1 - x2) or abs(y1 - y2) or abs(z1 - z2) > distance:
        return False
    return (x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2 > distance^2 # sqrt is slow

3D Euclidean distance is $d = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2+(z_1 - z_2)^2}$ but square roots are slow. Since we're just comparing relative distance, I think we could drop the sqrt and consider "distance" to be sums of squares.

The minimum calculation time would be if no 2-point matches are found. If there are 2000 points in the target, this would be 2000 * 2000 distance calculations, though many would be disqualified by a subtraction, and the results of previous calculations could be stored so you only have to do $2000 \choose 2$ = 1,999,000.

Actually, since you will need to calculate all of these anyway, whether you find matches or not, and since you only care about nearest neighbors for this step, if you have the memory it's probably better to pre-compute these values using an optimized algorithm. Something like a Delaunay or Pitteway triangulation, where every point in the target is connected to its nearest neighbors. Store those in a table, then look them up for each point when trying to fit the source triangle to one of the target triangles.

There's a lot of calculations involved, but it should be relatively fast since it's only operating on the data, which is sparse, instead of multiplying lots of meaningless zeros together like cross-correlation of volumetric data would involve. This same idea would work for the 2D case if you found the centers of the dots first and stored them as a set of coordinates.

Answer 2

I would add @mirror2image description on the alternative solution beside RANSAC, you may consider ICP algorithm (iterative closest point), a description can be found here!

I think the next challenge in using this ICP is to define your own cost function, and the starting pose of the target plane with respect to the 3d cloud point data. Some practical approach is to introduce some random noise in the data during the iteration to avoid converging to the false minima. This is the heuristic part I guess you need to design.

Update:

The steps in simplified form is:

1. Find the closest point for each input point.
2. Compute the transform from the input to the target, and then move the input points using the transform.
3. Compute the similarity function (e.g. the distance for each input point w.r.t to its corresponding pair target point).
4. Check the stopping condition.

Iterate the step 1-4.

There is available library that you can consider here! (I have not tried it yet), there is one section on the registration part (including other methods).