Tuomas Siipola Articles Projects

Aligning point patterns with Kabsch–Umeyama algorithm

  • Published
  • Last updated
  • Revisions

Kabsch–Umeyama algorithm is a method for aligning and comparing the similarity between two sets of points. It finds the optimal translation, rotation and scaling by minimizing the root-mean-square deviation (RMSD) of the point pairs. The algorithm for finding the optimal rotation was first described by Kabsch (1976, 1978). In this article we're taking a look at a similar method proposed by Umeyama (1991) which, in addition to rotation, also supports scaling.

Kabsch–Umeyama algorithm is efficient and easy to implement because it only requires a few matrix operations. It has many applications in different domains including physics simulation. However, the algorithm requires a predefined mapping between the sets of points. If you don't have this mapping, your data is noisy or incomplete, or you need a non-rigid transformation, methods for point set registration may work better for you.

Example

Before diving into the details, let's test out the algorithm in practise. Below we'll attempt to align the Little Dipper and Big Dipper constellations:

RMSD: 134.14
Step 1/5
JavaScript required

Algorithm

Inputs of the algorithm are two paired sets of \(n\) points represented as \(m\)-dimensional vectors. Let's say \(A = (\textbf{a}_1, \dots, \textbf{a}_n)\) are the reference points which we want to align points \(B = (\textbf{b}_1, \dots, \textbf{b}_n)\) with. We can write out the points as matrices which looks like the following in 2-dimensional case:

\[ A = \begin{bmatrix} x_1 & y_1 \\ x_2 & y_2 \\ \vdots & \vdots \\ x_n & y_n \end{bmatrix} \quad B = \begin{bmatrix} i_1 & j_1 \\ i_2 & j_2 \\ \vdots & \vdots \\ i_n & j_n \end{bmatrix} \]

The algorithm starts by finding the centroids of the points which are their arithmetic means in this case:

\[ \pmb\mu_A = \frac{1}{n}\sum_{i=1}^n \textbf a_i \quad \pmb\mu_B = \frac{1}{n}\sum_{i=1}^n \textbf b_i \]

We also need to the following variance to calculate scaling:

\[\sigma_A^2 = \frac{1}{n}\sum_{i=1}^n ||\textbf a_i - \pmb\mu_A||^2\]

Next, we calculate the following covariance matrix:

\[H = \frac{1}{n}\sum_{i=1}^n (\textbf{a}_i - \pmb\mu_A)^T(\textbf{b}_i - \pmb\mu_B)\]

and compute its singular value decomposition \(H = UDV^T\).

Reflection is sometimes allowed in addition to translation, rotation and scaling. In order to detect and prevent reflection we need the following matrix:

\[ \begin{aligned} d &= \text{sgn}(\text{det}(U)\,\text{det}(V^T)) \\ S &= \text{diag}(\underbrace{1, \dots, 1}_{m-1}, d) \end{aligned} \]

Finally, we can calculate the optimal rotation matrix \(R\) and scale factor \(c\):

\[ \begin{aligned} R &= USV^T \\ c &= \frac{\sigma_A^2}{\text{tr}(DS)} \\ \end{aligned} \]

With these at hand, we can calculate the aligned points \(B' = (\textbf{b}'_1, \dots, \textbf{b}'_n)\) where:

\[ \textbf{b}'_i = \pmb\mu_A + cR(\textbf{b}_i - \pmb\mu_B) \]

Sometimes it's nicer to have a single translation vector \(\textbf{t}\):

\[\textbf{t} = \pmb\mu_A - cR\pmb\mu_B\]

which can then be used to define the alignment as:

\[\textbf{b}'_i = \textbf{t} + cR\textbf{b}_i\]

Implementation

Translating the previous steps into code should be straightforward with linear algebra libraries or languages with such things built-in. Here's an implementation in Python with NumPy:

def kabsch_umeyama(A, B):
    assert A.shape == B.shape
    n, m = A.shape

    EA = np.mean(A, axis=0)
    EB = np.mean(B, axis=0)
    VarA = np.mean(np.linalg.norm(A - EA, axis=1) ** 2)

    H = ((A - EA).T @ (B - EB)) / n
    U, D, VT = np.linalg.svd(H)
    d = np.sign(np.linalg.det(U) * np.linalg.det(VT))
    S = np.diag([1] * (m - 1) + [d])

    R = U @ S @ VT
    c = VarA / np.trace(np.diag(D) @ S)
    t = EA - c * R @ EB

    return R, c, t

Here's how to align the constellations as shown above:

A = np.array([[ 23, 178],
              [ 66, 173],
              [ 88, 187],
              [119, 202],
              [122, 229],
              [170, 232],
              [179, 199]])
B = np.array([[232, 38],
              [208, 32],
              [181, 31],
              [155, 45],
              [142, 33],
              [121, 59],
              [139, 69]])

R, c, t = kabsch_umeyama(A, B)
# R = [[-0.81034281,  0.58595608]
#      [-0.58595608, -0.81034281]]
# c = 1.46166131
# t = [271.3345951, 396.07800317]

B = np.array([t + c * R @ b for b in B])
# B = [[ 29.08878779, 152.36814188]
#      [ 52.37669337, 180.03008629]
#      [ 83.50028582, 204.33920503]
#      [126.28647155, 210.02515345]
#      [131.40664707, 235.37261559]
#      [178.54823113, 222.56285654]
#      [165.79288328, 195.30194121]]

References