Categories &

Functions List

Function Reference: procrustes

statistics: d = procrustes (X, Y)
statistics: d = procrustes (X, Y, param1, value1, …)
statistics: [d, Z] = procrustes (…)
statistics: [d, Z, transform] = procrustes (…)

Procrustes Analysis.

d = procrustes (X, Y) computes a linear transformation of the points in the matrix Y to best conform them to the points in the matrix X by minimizing the sum of squared errors, as the goodness of fit criterion, which is returned in d as a dissimilarity measure. d is standardized by a measure of the scale of X, given by

  • sum (sum ((X - repmat (mean (X, 1), size (X, 1), 1)) .^ 2, 1))

i.e., the sum of squared elements of a centered version of X. However, if X comprises repetitions of the same point, the sum of squared errors is not standardized.

X and Y must have the same number of points (rows) and procrustes matches the i-th point in Y to the i-th point in X. Points in Y can have smaller dimensions (columns) than those in X, but not the opposite. Missing dimensions in Y are added with padding columns of zeros as necessary to match the the dimensions in X.

[d, Z] = procrustes (X, Y) also returns the transformed values in Y.

[d, Z, transform] = procrustes (X, Y) also returns the transformation that maps Y to Z.

transform is a structure with fields:

cthe translation component
Tthe orthogonal rotation and reflection component
bthe scale component

So that Z = transform.b * Y * transform.T + transform.c

procrustes can take two optional parameters as Name-Value pairs.

[…] = procrustes (…, 'Scaling', false) computes a transformation that does not include scaling, that is transform.b = 1. Setting 'Scaling' to true includes a scaling component, which is the default.

[…] = procrustes (…, 'Reflection', false) computes a transformation that does not include a reflection component, that is transform.T = 1. Setting 'Reflection' to true forces the solution to include a reflection component in the computed transformation, that is transform.T = -1.

[…] = procrustes (…, 'Reflection', 'best') computes the best fit procrustes solution, which may or may not include a reflection component, which is the default.

See also: cmdscale

Source Code: procrustes

Example: 1

Create some random points in two dimensions

 n = 10;
 randn ('seed', 1);
 X = normrnd (0, 1, [n, 2]);

Those same points, rotated, scaled, translated, plus some noise

 S = [0.5, -sqrt(3)/2; sqrt(3)/2, 0.5]; # rotate 60 degrees
 Y = normrnd (0.5*X*S + 2, 0.05, n, 2);

Conform Y to X, plot original X and Y, and transformed Y

 [d, Z] = procrustes (X, Y);
 plot (X(:,1), X(:,2), 'rx', Y(:,1), Y(:,2), 'b.', Z(:,1), Z(:,2), 'bx');
plotted figure

Example: 2

Find Procrustes distance and plot superimposed shape

 X = [40 88; 51 88; 35 78; 36 75; 39 72; 44 71; 48 71; 52 74; 55 77];
 Y = [36 43; 48 42; 31 26; 33 28; 37 30; 40 31; 45 30; 48 28; 51 24];
 plot (X(:,1),X(:,2),'x');
 hold on
 plot (Y(:,1),Y(:,2),'o');
 xlim ([0 100]);
 ylim ([0 100]);
 legend ('Target shape (X)', 'Source shape (Y)');
 [d, Z] = procrustes (X, Y)
d = 0.2026
Z =

   39.769   87.509
   50.562   86.801
   35.549   72.163
   37.313   73.991
   40.873   75.850
   43.552   76.796
   48.058   75.977
   50.783   74.229
   53.541   70.684
 plot (Z(:,1), Z(:,2), 's');
 legend ('Target shape (X)', 'Source shape (Y)', 'Transformed shape (Z)');
 hold off
plotted figure

Example: 3

Apply Procrustes transformation to larger set of points

Create matrices with landmark points for two triangles

 X = [5, 0; 5, 5; 8, 5];   # target
 Y = [0, 0; 1, 0; 1, 1];   # source

Create a matrix with more points on the source triangle

 Y_mp = [linspace(Y(1,1),Y(2,1),10)', linspace(Y(1,2),Y(2,2),10)'; ...
         linspace(Y(2,1),Y(3,1),10)', linspace(Y(2,2),Y(3,2),10)'; ...
         linspace(Y(3,1),Y(1,1),10)', linspace(Y(3,2),Y(1,2),10)'];

Plot both shapes, including the larger set of points for the source shape

 plot ([X(:,1); X(1,1)], [X(:,2); X(1,2)], 'bx-');
 hold on
 plot ([Y(:,1); Y(1,1)], [Y(:,2); Y(1,2)], 'ro-', 'MarkerFaceColor', 'r');
 plot (Y_mp(:,1), Y_mp(:,2), 'ro');
 xlim ([-1 10]);
 ylim ([-1 6]);
 legend ('Target shape (X)', 'Source shape (Y)', ...
         'More points on Y', 'Location', 'northwest');
 hold off

Obtain the Procrustes transformation

 [d, Z, transform] = procrustes (X, Y)
d = 0.044118
Z =

   5.0000   0.5000
   4.5000   4.5000
   8.5000   5.0000

transform =

  scalar structure containing the fields:

    T =

      -0.1240   0.9923
       0.9923   0.1240

    b = 4.0311
    c =

       5.0000   0.5000
       5.0000   0.5000
       5.0000   0.5000

Use the Procrustes transformation to superimpose the more points (Y_mp) on the source shape onto the target shape, and then visualize the results.

 Z_mp = transform.b * Y_mp * transform.T + transform.c(1,:);
plotted figure 

 figure
 plot ([X(:,1); X(1,1)], [X(:,2); X(1,2)], 'bx-');
 hold on
 plot ([Y(:,1); Y(1,1)], [Y(:,2); Y(1,2)], 'ro-', 'MarkerFaceColor', 'r');
 plot (Y_mp(:,1), Y_mp(:,2), 'ro');
 xlim ([-1 10]);
 ylim ([-1 6]);
 plot ([Z(:,1); Z(1,1)],[Z(:,2); Z(1,2)],'ks-','MarkerFaceColor','k');
 plot (Z_mp(:,1),Z_mp(:,2),'ks');
 legend ('Target shape (X)', 'Source shape (Y)', ...
         'More points on Y', 'Transformed source shape (Z)', ...
         'Transformed additional points', 'Location', 'northwest');
 hold off
plotted figure

Example: 4

Compare shapes without reflection

 T = [33, 93; 33, 87; 33, 80; 31, 72; 32, 65; 32, 58; 30, 72; ...
      28, 72; 25, 69; 22, 64; 23, 59; 26, 57; 30, 57];
 S = [48, 83; 48, 77; 48, 70; 48, 65; 49, 59; 49, 56; 50, 66; ...
      52, 66; 56, 65; 58, 61; 57, 57; 54, 56; 51, 55];
 plot (T(:,1), T(:,2), 'x-');
 hold on
 plot (S(:,1), S(:,2), 'o-');
 legend ('Target shape (d)', 'Source shape (b)');
 hold off
 d_false = procrustes (T, S, 'reflection', false);
 printf ("Procrustes distance without reflection: %f\n", d_false);
Procrustes distance without reflection: 0.342463
 d_true = procrustes (T, S, 'reflection', true);
 printf ("Procrustes distance with reflection: %f\n", d_true);
Procrustes distance with reflection: 0.020428
 d_best = procrustes (T, S, 'reflection', 'best');
 printf ("Procrustes distance with best fit: %f\n", d_true);
Procrustes distance with best fit: 0.020428
plotted figure