function procrustes
close all
Shape analysis on a beach holiday
Procrustes analysis is a method for analysing sets of shapes (see [1]). Let's suppose we pick up a pebble from the beach and want to know how closely its shape matches the outline of a frisbee. Here is a plot of the frisbee and the pebble.
figure(1)
t = chebfun('x',[0,2*pi]);
f = 3*(1.5*cos(t) + 1i*sin(t)); %frisbee
g = exp(1i*pi/3)*(1+cos(t)+1.5i*sin(t)+ .125*(1+1.5i)*sin(3*t).^2); %pebble
plot(f,'r','LineWidth',2), hold on, axis equal, plot(g,'k','LineWidth',2)
title('Frisbee and pebble','FontSize',14), hold off
Two shapes are equivalent if one can be obtained from the other by translating, scaling and rotating. Before comparison we thus:
-
Translate the shapes so they have mean zero.
-
Scale so the shapes have Root Mean Squared Distance (RMSD) to the origin of $1$.
-
Rotate to align major axis.
Here is how the frisbee and the pebble compare after each stage.
function [f,g]=ShapeAnalysis(f,g)
% SHAPEANALYSIS(F,G) Plots the parameterised curves before and after
% each stage of translating, scaling and aligning. Outputs are
% parameterised curves ready for Procrustes shape analysis.
LW = 'LineWidth'; FS = 'FontSize';
% Plot original
subplot(2,2,1)
plot(f,'r',LW,2), hold on, axis equal, plot(g,'k',LW,2)
title('Original',FS,14)
% Translate mean to 0.
f = f-mean(f); g = g-mean(g);
subplot(2,2,2)
plot(f,'r',LW,2), hold on, axis equal, plot(g,'k',LW,2)
title('After translation',FS,14)
% Scale so RMSD is 1.
f = f/norm(f); g = g/norm(g);
subplot(2,2,3)
plot(f,'r',LW,2), hold on, axis equal, plot(g,'k',LW,2)
title('After scaling',FS,14)
% Align major axis.
subplot(2,2,4)
% Find argument of major axis.
[ignored,fxmax] = max(abs(f));
[ignored,gxmax] = max(abs(g));
rotf = angle(f(fxmax)); rotg = angle(g(gxmax));
% Rotate both so major axis lies on the positive real axis.
x = chebfun('x',[0,2*pi]);
f = exp(-1i*rotf)*f(mod(x+fxmax,2*pi));
g = exp(-1i*rotg)*g(mod(x+gxmax,2*pi));
plot(f,'r',LW,2), hold on, axis equal, plot(g,'k',LW,2)
title('After aligning',FS,14), hold off
end
[f,g] = ShapeAnalysis(f,g);
To calculate the Procrustes distance in the discrete case we would measure the error between the two shapes at a finite number of reference points and compute the vector $2$-norm. In Chebfun we calculate the continuous analogue:
norm(f-g)
ans = 0.072347575424997
A little warning
In the discrete version of Procrustes shape analysis statisticians choose reference points on the two shapes (to compare). They then work out the difference between corresponding reference points. The error computed depends on this correspondence. A different correspondence gives a different error. In the continuous case this correspondence becomes the parameterisation. A different parameterisation of the two curves gives a different error. This continuous version of Procrustes (as implemented in this example) is therefore more of an 'eyeball' check than a robust statistical analysis.
A shape and its reflection
At the beach shapes reflect on the surface of the sea. An interesting question is: How close, in shape, is a pebble to its reflection? Here is a plot of a pebble and its reflection.
figure(2)
% pebble
f = exp(1i*pi/3)*(1+cos(t)+1.5i*sin(t)+.125*(1+1.5i)*sin(3*t).^2);
% reflection
g = exp(-1i*pi/3)*(1+cos(2*pi-t)-1.5i*sin(2*pi-t)+.125*(1-1.5i)*sin(3*(2*pi-t)).^2);
plot(f,'r','LineWidth',2), hold on, axis equal, plot(g,'k','LineWidth',2)
title('Pebble and its reflection','FontSize',14), hold off
Here is how the pebble and its reflection compare after each stage of translating, scaling and rotating.
[f,g] = ShapeAnalysis(f,g);
Now we calculate the continuous Procrustes distance.
norm(f-g)
ans = 0.097593759012228
Comparing this result to the Procrustes distance of the pebble and a frisbee shows that the pebble is closer in shape to a frisbee than to its own reflection!
end