Orthogonal Procrustes Analysis with isotropic scaling

Orthogonal Procrustes Analysis (OPA) refer to the problem that relates two datasets, with known correspondence, through a transformation. In the context of computer-assisted intervention, this problem is often encountered as either a registration or calibration between 2 sets of 3D points. It is also known as landmark-based registration, paired fiducial registration, among others.

For example, $x_i$ may refer to $i^{th}$ fiducial landmark in a CT volume and $y_i$ be the corresponding (refer to the same anatomical feature) fiducial in a MRI volume. Given sufficient number of these corresponding landmarks, these two 3D volumes can be aligned/registered.

The transformation one seeks typically comprises of a scaling factor, followed by rotation and translation. Formally, we want to minimize the Fiducial Registration Error (FRE) between the point set $Y$ and $X$ after registration:

Many solutions exists since 1940’, but in the engineering field, the most commonly cited solution is either the quaternion based [1] or the singular-value-decomposition (SVD) based solution [2]. While the quaternion based solution only works in 3D, it is not prone to erroneous reflection (mirroring). SVD-based solution can be generalized to mD, but requires an extra step to prevent reflection in rotation. The following solution is based on [2] for solving the optimal rotation and translational component of the transformation, but uses the formula from [1] for solving the isotropic scaling factor. Note that points are stored as column vectors.

function [R, t, a] = OPA( X, Y )
% Solves Y = R a X + t, least-square solution
% INPUTS:  X : (mxn) mD coordinates,
%          Y : (mxn) mD coordinates
% OUTPUTS: R : mxm orthonormal rotation
%          t : mx1 translation
%          a : a scalar scaling factor
%
% Elvis C.S. Chen
% Aug. 24, 2020
[m,n] = size(X); e = ones(n,1); II = eye(n)-((e*e')./n);
[U,~,V] = svd(Y*II*X'); D = eye(m); D(m,m)=det(U*V');
R = U * D * V';                                      % rotation
a = sqrt( sum(sum((Y*II).^2))/sum(sum((X*II).^2)) ); % scaling
t = mean( Y-R*a*X, 2);                               % translation

The code above is optimize for length, not necessary for computational efficiency.

References

[1] Horn, B. K. P. (1987). Closed-form solution of absolute orientation using unit quaternions. Journal of the Optical Society of America A, 4(4), 629–642. https://doi.org/10.1364/JOSAA.4.000629

[2] Arun, K. S., Huang, T. S., & Blostein, S. D. (1987). Least-Squares Fitting of Two 3-D Point Sets. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-9(5), 698–700. https://doi.org/10.1109/TPAMI.1987.4767965