Superposition

3-D Rotation Problem: Given S = s₀, s₁, ..., s_n-1, and T = t₀, t₁, ... t_n-1, where each s_i and each t_i 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 - t_i))².

 

Note, it is given that s_i ~ t_i,  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).

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[0]-ti[0], xp = si[0]+ti[0],
   ym = si[1]-ti[1], yp = si[1]+ti[1],
   zm = si[2]-ti[2], zp = si[2]+ti[2];

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[0], vec = val_vec[1];
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].