### Superposition

 LA home Computing Algorithms  glossary  Numerical   Num.Errors   Polynomials   Stirling   Mean&S.D.   Segmentation   Superposition   Integration   Matrices   Eigen v.   Superposition
3-D Rotation Problem: Given S = s0, s1, ..., sn-1, and T = t0, t1, ... tn-1, where each si and each ti is a 3-D vector, ⟨?x, ?y, ?z⟩, find a rotation, R, by some angle θ, about some axis, u, through the origin, ⟨0, 0, 0⟩, so as to minimise the sum of squared errors between S' = map R S and T, that is  i (norm(s'i - ti))2.

Note, it is given that si ~ ti 0 ≤ i < n.
To find the general rigid-body transformation of S that gives the least-squares superposition of S on T (a version of the Total Least-Squares (TLS) problem), first make the centre of gravity (CoG) of T the origin, then translate S's CoG (carrying S with it) to that origin, and finally solve the rotation problem above.

This problem appears in the superposition of two protein structures in Bioinformatics.

Kearsley's Method: This uses the quaternion representation of 3-D rotations, and an Eigen-value finding algorithm, such as the Jacobi algorithm, say. A pure rotation is represented by a constrained quaternion -- a general quaternion having four degrees of freedom and a 3-D rotation only three -- so a Lagrange multiplier also gets a look in. All this leads to an Eigen-value problem on a 4×4 matrix. The smallest Eigen-value, λ, is the least-square error, and the corresponding Eigen-vector yields the required rotation (quaternion).
 S: 0 0 0 1 0 0 0 1 0 0 0 1 T: 0 0 0 1 0 0.1 0 0 -1 0 1 0
 -- L.A., Aug. 2011.
output
```
xm2 = ∑i xm*xm; xp2 = ∑i xp*xp;
ym2 = ∑i ym*ym; yp2 = ∑i yp*yp;
zm2 = ∑i zm*zm; zp2 = ∑i zp*zp;
xmym = ∑i xm*ym;
xmyp = ∑i xm*yp;
xmzm = ∑i xm*zm;
xmzp = ∑i xm*zp;
xpym = ∑i xp*ym;
xpyp = ∑i xp*yp;
xpzm = ∑i xp*zm;
xpzp = ∑i xp*zp;
ymzm = ∑i ym*zm;
ymzp = ∑i ym*zp;
ypzm = ∑i yp*zm;
ypzp = ∑i yp*zp;
where
xm = si-ti, xp = si+ti,
ym = si-ti, yp = si+ti,
zm = si-ti, zp = si+ti;

M = [ [xm2+ym2+zm2,   ypzm-ymzp,   xmzp-xpzm,   xpym-xmyp],
[  ypzm-ymzp, xm2+yp2+zp2,   xmym-xpyp,   xmzm-xpzp],
[  xmzp-xpzm,   xmym-xpyp, xp2+ym2+zp2,   ymzm-ypzp],
[  xpym-xmyp,   xmzm-xpzp,   ymzm-ypzp, xp2+yp2+zm2] ];

val_vec = Jacobi(M);  // Eigen-values and -vectors
val = val_vec, vec = val_vec;
small = 0;
for(i = 1; i < val.length; i ++)
if(val[i] < val[small]) small=i;
return rotationQ( vec[small] );
```
Thanks to Arun K for bringing this method to my attention.
Search for [Kearsley superposition] in the [Bib].