The Kabsch algorithm, also known as the Kabsch-Umeyama algorithm, named after Wolfgang Kabsch and Shinji Umeyama, is a method for calculating the optimal rotation matrix that minimizes the RMSD (root mean squared deviation) between two paired sets of points. It is useful for point-set registration in computer graphics, and in cheminformatics and bioinformatics to compare molecular and protein structures (in particular, see root-mean-square deviation (bioinformatics)). The algorithm only computes the rotation matrix, but it also requires the computation of a translation vector. When both the translation and rotation are actually performed, the algorithm is sometimes called partial Procrustes superimposition (see also orthogonal Procrustes problem).
Description Let P and Q be two sets, each containing N points in R n {\displaystyle \mathbb {R} ^{n}} . We want to find the transformation from Q to P. For simplicity, we will consider the three-dimensional case ( n = 3 {\displaystyle n=3} ). The sets P and Q can each be represented by N × 3 matrices with the first row containing the coordinates of the first point, the second row containing the coordinates of the second point, and so on, as shown in this matrix:
( x 1 y 1 z 1 x 2 y 2 z 2 ⋮ ⋮ ⋮ x N y N z N ) {\displaystyle {\begin{pmatrix}x_{1}&y_{1}&z_{1}\\x_{2}&y_{2}&z_{2}\\\vdots &\vdots &\vdots \\x_{N}&y_{N}&z_{N}\end{pmatrix}}}
The algorithm works in three steps: a translation, the computation of a covariance matrix, and the computation of the optimal rotation matrix.
Translation Both sets of coordinates must be translated first, so that their centroid coincides with the origin of the coordinate system. This is done by subtracting the centroid coordinates from the point coordinates.
Computation of the covariance matrix The second step consists of calculating a matrix H. In matrix notation,
H = P T Q {\displaystyle H=P^{\mathsf {T}}Q\,}
or, using summation notation,
H i j = ∑ k = 1 N P k i Q k j , {\displaystyle H_{ij}=\sum _{k=1}^{N}P_{ki}Q_{kj},}
which is a cross-covariance matrix when P and Q are seen as data matrices.
Computation of the optimal rotation matrix It is possible to calculate the optimal rotation R based on the matrix formula
R = ( H T H ) 1 2 H − 1 , {\displaystyle R=\left(H^{\mathsf {T}}H\right)^{\frac {1}{2}}H^{-1},}
but implementing a numerical solution to this formula becomes complicated when all special cases are accounted for (for example, the case of H not having an inverse). If singular value decomposition (SVD) routines are available the optimal rotation, R, can be calculated using the following algorithm. First, calculate the SVD of the covariance matrix H,
H = U Σ V T {\displaystyle H=U\Sigma V^{\mathsf {T}}}
where U and V are orthogonal and Σ {\displaystyle \Sigma } is diagonal. Next, record if the orthogonal matrices contain a reflection,
d = det ( U V T ) = det ( U ) det ( V ) . {\displaystyle d=\det \left(UV^{\mathsf {T}}\right)=\det(U)\det(V).}
… excerpt ends here. Continue reading the full article.
