4

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!

froost
  • 41
  • 1
  • 4
  • This is a very difficult problem. With differing numbers of points it is entirely possible that you can find transformations such that the set with less points is a subset of the other. Furthermore there could be many such transformations depending on your data and then how do you decide which is best? Is something simple suitable such as overlaying the centroid of each set and then finding the optimal rotation that minimised RMS ?? Or is more accuracy needed ? – mathematician1975 Jun 23 '12 at 08:41
  • @mathematician1975 : TBH, the problem is not as open ended as it sounds. The number of points are usually similar, the larger set might have max of twice the points, and my goal is to find "a" (not optimal) superposition that aligns some >30-50% of the points in a way they overlap. If a perfect rotation exists (the smaller set is a subset after transformation), that is the perfect case I'm looking for. Otherwise, it should decrease the size until an acceptable match is made. Centroids might be OK, but think: a sphere & a hemisphere. Their centroids are different in the perfect alignment. – froost Jun 23 '12 at 10:36
  • 1
    There isn't a clear migration path. This is on topic here, but might be too academic to get the traction you need. You might want to check on the Metas for some of the relevant migration paths. Contacted a [cstheory.se] mod and he was of the opinion that it wouldn't be a good fit there. Might get more traction on [gamedev.se] or [cs.se], but I'm not sure. –  Jun 24 '12 at 16:52
  • are the point sets arbitrary? or are they samples of some underlying geometric structure that can be exploited by the solver? – atb Jun 25 '12 at 00:17
  • the points are from protein structures, so there is an underlying structure, however I think not something that can easily be exploited. – froost Jun 27 '12 at 14:53

3 Answers3

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.

datjko
  • 383
  • 1
  • 10
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 · xi = bi
  3. Extend to multiple vectors: A · X = B
  4. Transpose: (A · X)T = BT
  5. Simplify: XT · AT = BT
  6. Substitute P = XT, Q = AT and R = BT. The result is: P · Q = R
  7. Apply the formula for LLS: Q ≈ (PT · P)-1 · PT · R.
  8. Substitute back: AT ≈ (X · XT)-1 · X · BT
  9. Solve for A, and simplify: A ≈ B · XT · (X · XT)-1

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

  • B · XT = ∑bi·xiT
  • X · XT = ∑xi·xiT
  • A ≈ (∑bi·xiT) · (∑xi·xiT)-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.

Markus Jarderot
  • 86,735
  • 21
  • 136
  • 138
0

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

Community
  • 1
  • 1
Tolga Birdal
  • 491
  • 6
  • 15