# Procrustes analysis

Standard

Procrustes Analysis is a method to compare shapes of two or more objects, by first performing superimposition, i.e. rotating, translating, and uniformly scaling the objects so that they fit each other best. Here’s an example of it, and sample C++ OpenCV code for it.

Basically we have a set of points X (green dots in the picture), and another set of points Y (blue ones). Points in Y are actually points in X that have undergone rotation, scaling, and translation transformations, plus a bit of jittering to randomize the final result.To perform procrustes analysis on X and Y, we follow the following steps:

1. Recenter all the points and normalise them. Basically subtract each points with their means, and perform normalisation on them
2. Perform SVD to obtain the U, S, V components. Since U S V components can be viewed as rotation x scaling x another rotation, we could derive or rotation and uniform scale transformation from it.
3. Since we now have more or less all the transformations on our hand, we could find Y’ by subjecting Y (or normalised Y) to those transformations.

Here’s the result of applying Procrustes method on it. The red dots are Y’, which is Y after applying the transformations found from Procrustes analysis.

Here’s a snippet of the code of the Procrustes function that I ported from http://stackoverflow.com/questions/18925181/procrustes-analysis-with-numpy

```float Procrustes::procrustes( const Mat& X, const Mat& Y ){
/* Recenter the points based on their mean ... */
Scalar mu_x = cv::mean(X);
Mat X0      = X - Mat(X.size(), X.type(), mu_x);

Scalar mu_y = cv::mean(Y);
Mat Y0      = Y - Mat(Y.size(), Y.type(), mu_y);

/* ... and normalize them */
float ss_X      = sumSquared( X0 );
float norm_X    = sqrt( ss_X );
X0              /= norm_X;

float ss_Y      = sumSquared( Y0 );
float norm_Y    = sqrt( ss_Y );
Y0              /= norm_Y;

/* Pad with zeros is Y has less points than X */
if( Y.rows < X.rows )
vconcat( Y0, Mat::zeros( X.rows - Y.rows, 1, Y.type()), Y0 );

/* Perform SVD */
Mat A = X0.reshape(1).t() * Y0.reshape(1);
Mat U, s, Vt;
SVDecomp( A, s, U, Vt );

Mat V           = Vt.t();
this->rotation  = V * U.t();

if( !bestReflection ) {
bool have_reflection = determinant( this->rotation ) < 0;
if( bestReflection != have_reflection ) {
V.colRange( V.cols-1, V.cols) *= -1;
s.rowRange( s.rows-1, s.rows) *= -1;
this->rotation = V * U.t();
}
}

/* Rotate Y0 first */
Mat rotated_Y0;
cv::transform( Y0, rotated_Y0, rotation.t() );

/* Trace of eigenvalues is basically the scale */
float trace_TA = sum( s )[0];

if( scaling ) {
scale   = trace_TA * norm_X / norm_Y;
error   = 1 - trace_TA * trace_TA;
Yprime  = norm_X * trace_TA * rotated_Y0 + mu_x;
}
else {
error   = 1 + ss_Y / ss_X - 2 * trace_TA * norm_X / norm_Y;
Yprime  = norm_Y * rotated_Y0 + mu_x;
}

if( Y.rows < X.rows )
rotation = rotation.rowRange(0, Y.rows );

translation = Mat(1, 1, CV_32FC2, mu_x).reshape(1) - scale * Mat(1, 1, CV_32FC2, mu_y).reshape(1) * rotation;

return error;
}
```

You can get the sample code to generate the data and calculate the Procrustes analysis here: https://github.com/subokita/Sandbox/tree/master/Procrustes/Procrustes