Procrustes analysis is used for comparing shapes of objects, the word generalised being used when there are more than 2 objects involved. The shape of an object is defined ([5]) in a mathematical context as all the geometrical information that remains after location, scale and rotational effects are filtered out; i.e. an obect's shape is invariant under the Euclidean similarity transformations of translation, scaling and rotation. For analysis purposes shape is described by a finite number of points on each object's surface which are called landmarks. So a landmark is a point of correspondence on each object that matches between and within populations. Three basic types of landmarks are generally used :

(1) anatomical landmarks - points assigned by an "expert" which correspond between organisms in some biologically meaningful way e.g. the points on the skull in figure 1.1

(2) mathematical landmarks - points located on an object according to some mathematical or geometrical property of the object

(3) pseudo-landmarks - points located either around the outline (2D) or between anatomical or mathematical landmarks

The configuration is the set of landmarks on a particular object and the configuration matrix X is the kxm matrix of Cartesian coordinates of the k landmarks in m dimensions (m = 2 or 3). The configurations being considered will all be ordered sets, so that the ith row of configuration matrix X will correspond to the coordinates of the ith landmark on X, which will in turn correspond to the ith landmark on every other object being compared.

A rotation of a configuration matrix X_kxm is given by post-multiplication of X by a rotation matrix i.e. an mxm matrix R satisfying R^TR = I and |R| = 1. The set of all mxm rotation matrices is denoted by SO(m). A translation is obtained by adding a constant to each entry in X, while an isotropic scaling is obtained by multiplying X by a positive real number. The Euclidean similarity transformations of X are all the possible rotated, translated and isotropically scaled versions of X and so may be represented by the following set :

{aXR + 1_kv^T : a R⁺, RSO(m), v R^m}

where a R⁺ is the scale, R is an mxm rotation matrix, v is an mx1 translation vector (v₁, v₂, v₃)^T say (m=3), and 1_k is the kx1 vector of ones ( so 1_kv^T is the kxm matrix with each row equal to v^T).

A size measure g(X) of a configuration matrix X is a positive real valued function of X such that

g(aX) = ag(X) for any positive scalar a

The most commonly used size measure is the centroid size given by

C = I_k - (1_k1_k^T) is the centring matrix,

Looking closer at CX it is clear that this has the effect on X of translating X so that it is centred (at the origin) i.e. the sum of its rows (landmark coordinates) is zero :

Full generalised Procrustes analysis is a way of finding the similarity transforms which need to be applied to two or more configuration martices X_i so that they are as close to each other as possible in terms of the Eudlidean norm. It may be thought of in 3D (or 2D) as rotating, scaling and translating all the X_i so that the sum of the volumes (or areas) of the smallest spheres (or circles) which enclose each set of corresponding vertices is at a minimum; or more algebraically, as minimising the sum of the Euclidean distances/norms of each shape from every other shape. It is a generalisation of full Ordinary Procrustes Analysis (OPA) which finds the similarity transformations to be applied to one configuration matrix X₁ which minimise its Euclidean distance from a second configuration matrix X₂. i.e. it finds the similarity parameters b, R and v which minimise :

D²_OPA(X₁, X₂) =

where RSO(m), b > 0 is a scale parameter and v is a location (translation) vector.

The minimum of the above equation is denoted by OSS(X₁, X₂) which stands for Ordinary (Procrustes) Sum of Squares.

For centred configuration matrices X₁ and X₂, it can be shown that D²_OPA(X₁, X₂) = OSS(X₁, X₂) when v = 0, R = UV^T, where X₂X₁^T = VDU^T is the SVD of X₂X₁^T, and b = tr(X₂^TX₁R)/tr(X₁^TX₁).

So to perform a full OPA on two configuration matrices X₁^O and X₂^O , first centre the matrices to get X₁ and X₂, then find the SVD of X₂X₁^T (Java and Matlab functions exist to do this), then apply the rotation R and scale b given by the above result to the X₁ configuration matrix. The new configuration matrix for X₁ is called the full Procrustes fit (or full Procrustes coordinates) of X₁ onto X₂ and is denoted by X₁^P.Now suppose X₁, ..., X_n (n >1) are configuration matrices. Then the method of full generalised Procrustes analysis (GPA) finds the similarity parameters b_i, R_i and v_i which minimise (a quantity proportional to) the sum of squared norms of pairwise differences :

D²_GPA(X₁, ..., X_n) =

subject to a constraint on the size of the average,

S() = 1,

where R_iSO(m), b_i > 0, and S() is the centroid size of the average configuration

The minimum of D²_GPA(X₁, ..., X_n) subject to this constraint is denoted by G(X₁, ..., X_n) and is called the generalised (Procrustes) sum of squares. The new configuration matrices given by the minimising parameters are called the full Procrustes fits (or coordinates) of the X_i and are denoted by X_i^P.

For objects in more than 2 dimensions iterative procedures must be used for the matching process. For a practical implementation one may use the GPA algorithm of Gower(1975), modified by Ten Berge (1977) :

1. Translations. Centre the configuration matrices to remove location. Initially let

= , i = 1, ..., n.

2. Rotations. For the ith configuration let

.

Then the new is taken to be the ordinary Procrustes superimposition, involving only rotation, of the old on . The n figures are rotated in turn. This process is repeated until the Procrustes sum of squares D²_GPA(X₁, ..., X_n) (with b_i = 1, v_i = 0 i) cannot be reduced further.

3. Scaling. Let be the nxn correlation matrix of the vec(X_i^P) ( where vec(X_i^P) = vector made of concatenation of rows of X_i, and the usual roles of variable and observation labels are reversed) with eigenvector corresponding to the largest eigenvalue. Then (from Ten Berge(1977)) take

b_i = , repeated for all i.

4. Repeat steps 2 and 3 until the Procrustes sum of squares D²_GPA(X₁, ..., X_n) cannot be reduced further.