Algorithm for superimposition of 3d points

Question

I need to superimpose two groups of 3D points on top of each other; i.e. find rotation and translation matrices to minimize the RMSD (root mean square deviation) between their coordinates.

I currently use Kabsch algorithm, which is not very useful for many of the cases I need to deal with. Kabsch requires equal number of points in both data sets, plus, it needs to know which point is going to be aligned with which one beforehand. For my case, the number of points will be different, and I don't care which point corresponds to which in the final alignment, as long as the RMSD is minimized.

So, the algorithm will (presumably) find a 1-1 mapping between the subsets of two point sets such that AFTER rotation&translation, the RMSD is minimized.

I know some algorithms that deal with different number of points, however they all are protein-based, that is, they try to align the backbones together (some continuous segment is aligned with another continuous segment etc), which is not useful for points floating in space, without any connections. (OK, to be clear, some points are connected; but there are points without any connections which I don't want to ignore during superimposition.)

Only algorithm that I found is DIP-OVL, found in STRAP software module (open source). I tried the code, but the behaviour seems erratic; sometimes it finds good alignments, sometimes it can't align a set of few points with itself after a simple X translation.

Anyone know of an algorithm that deals with such limitations? I'll have at most ~10^2 to ~10^3 points if the performance is an issue.

---

To be honest, the objective function to use is not very clear. RMSD is defined as the RMS of the distance between the corresponding points. If I have two sets with 50 and 100 points, and the algorithm matches 1 or few points within the sets, the resulting RMSD between those few points will be zero, while the overall superposition may not be so great. RMSD between all pairs of points is not a better solution (I think).

Only thing I can think of is to find the closest point in set X for each point in set Y (so there will be exactly min(|X|,|Y|) matches, e.g. 50 in that case) and calculate RMSD from those matches. But the distance calculation and bipartite matching portion seems too computationally complex to call in a batch fashion. Any help in that area will help as well.

Thanks!

Answer 1

The best algorithm for discovering alignments through superimposition is Procrustes Analysis or Horn's method. Please follow this Stackoverflow link.

Answer 2

If you know which pairs of points correspond to each other, you can recover the transformation matrix with Linear Least Squares (LLS).

When considering LLS, you normally would want to find an approximation of x in A*x = b. With a transpose, you can solve for A instead of x.

1. Extend each source and target vector with "1", so they look like <x, y z, 1>
2. Equation: A · x_i = b_i
3. Extend to multiple vectors: A · X = B
4. Transpose: (A · X)^T = B^T
5. Simplify: X^T · A^T = B^T
6. Substitute P = X^T, Q = A^T and R = B^T. The result is: P · Q = R
7. Apply the formula for LLS: Q ≈ (P^T · P)^-1 · P^T · R.
8. Substitute back: A^T ≈ (X · X^T)^-1 · X · B^T
9. Solve for A, and simplify: A ≈ B · X^T · (X · X^T)^-1

(B · X^T) and (X · X^T) can be computed iteratively by summing up the outer products of the individual vector pairs.

• B · X^T = ∑b_i·x_i^T
• X · X^T = ∑x_i·x_i^T
• A ≈ (∑b_i·x_i^T) · (∑x_i·x_i^T)^-1

No matrix will be bigger than 4×4, so the algorithm does not use any excessive memory.

The result is not necessarily affine, but probably close. With some further processing, you can make it affine.

Answer 3

What you said looks like a "cloud to cloud registration" task. Take a look into http://en.wikipedia.org/wiki/Iterative_closest_point and http://www.willowgarage.com/blog/2011/04/10/modular-components-point-cloud-registration for example. You can play with your data in open source Point Cloud Library to see if it works for you.