# Aligning Point Patterns with Kabsch–Umeyama Algorithm

- Published
- 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:

## 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

- Kabsch, W. (1976). "A solution for the best rotation to relate two sets of vectors". Acta Crystallographica. A32 (5): 922–923. doi:10.1107/S0567739476001873.
- Kabsch, W. (1978). "A discussion of the solution for the best rotation to relate two sets of vectors". Acta Crystallographica. A34 (5): 827–828. doi:10.1107/S0567739478001680
- Umeyama, S. (1991). "Least-squares estimation of transformation parameters between two point patterns". IEEE Transactions on Pattern Analysis and Machine Intelligence. 13 (4): 376–380. doi:10.1109/34.88573.