Morphometrics, the measurement (metron) of shape (morphe), is a subfield of statistics with a history going back to the very beginnings of this discipline. For example, in 1888 Frances Galton introduced the correlation coefficient and applied it to a variety of morphological measurements on humans. In 1907 he invented a method to quantify facial shape that has later been termed two-point shape coordinates or Bookstein-shape coordinates (see below). The application of multivariate statistical techniques, which were basically developed in the first half of the 20th century, led to so-called multivariate morphometrics. But in the 1980s, morphometrics experienced a major revolution through the invention of coordinate-based methods, the discovery of the statistical theory of shape, and the computational realization of deformation grids (for historical reviews see Bookstein 1998; Rohlf and Marcus 1993; O’Higgins 2000; Adams et al. 2004; Slice 2005). The ubiquitous application of fast personal computers has ushered in a new era of data analysis, permitting the exploration and visualization of large high-dimensional data sets along with exact statistical tests based on resampling procedures.

This new morphometric approach has been termed geometric morphometrics as it preserves the geometry of the landmark configurations throughout the analysis and thus permits to represent statistical results as actual shapes or forms. Among several geometric approaches to morphometrics, the Procrustes method is the most widespread and best understood in its mathematical and statistical properties (Bookstein 1996; Small 1996; Dryden and Mardia 1998). Other frequently used morphometric methods are Euclidian distance matrix analysis (Lele and Richtsmeier 1991, 2001), elliptic Fourier analysis (Lestrel 1982), and non-label based approaches like voxel-based morphometry (e.g., Ashburner and Friston 2000), which is mainly applied in brain imaging. For more traditional morphometric approaches see Blackith and Reyment (1971), Marcus (1990), and Oxnard (1983), or the variety of methods applied in histology (Baak and Oort 1983) and stereology (Weibel 1979; Baddeley and Vedel Jensen 2004). For a rigorous statistical comparison of several of these methods see Rohlf (2000a, 2000b, 2003). This article focuses on Procrustes methods and the associated statistical and graphical toolkit (Table 1 provides definitions of some frequently used terms).

Table 1 Some key terms in geometric morphometrics and their synonyms

Geometric morphometrics is based on landmark coordinates. Bookstein (1991, p. 2) defines landmarks as loci that have names (‘bridge of the nose’, ‘tip of the chin’) as well as Cartesian coordinates. The names are intended to imply correspondence (biological homology) among forms. That is, landmark points not only have their own locations but also have the “same” locations in every other form of the sample and in the average of all the forms.

These coordinate data can come from a vast array of sources and can either be two- or three-dimensional. Two-dimensional coordinates are usually captured using a digitizing tablet or by measuring an image on the computer. Three-dimensional data can be captured directly using a coordinate digitizer such as a Microscribe or Polhemus, or may be measured on surface scans or volumetric scans. Volumetric data are based on image-slices from computed tomographic (CT) or magnet resonance imaging (MRI) scanners—or their high-resolution versions μCT and μMRI (Sensen and Hallgrimsson 2009). These slices contain gray-values that correspond to tissue densities and are concatenated to obtain a three-dimensional representation of an object. Surface scanners provide high-resolution 3D representations of an object’s surface using either laser or more traditional optical technology and may also include texture information. Surfaces can also be extracted from CT or MRI data. Most software packages allow landmark coordinates to be measured directly on these virtual surfaces or volumetric objects.

Traditional (i.e., non-geometric) morphometric approaches typically apply statistical techniques to a wide range of measurements, such as distances and distance ratios, angles, areas, and volumes. It is of course possible to compute the traditional interlandmark distances and angles from the landmark coordinates, but the original geometric relationship among the points may not be reconstructable from a sample of selected distance measurements. The raw landmark coordinates, however, cannot be subject to statistical analysis before separating shape information from overall size and “nuisance parameters” like position and orientation of the specimens in the digitizing space that are not relevant for the analysis.

The Geometry of Shape and Form

In morphometrics, the term shape is used to denote the geometric properties of an object that are independent of the object’s overall size, position, and orientation, whereas the form of an object comprises both its shape and size. Take, for example, the three landmarks A, B, and C in Fig. 1, constituting a triangle as the simplest two-dimensional configuration. The three interlandmark distances sufficiently describe the triangle’s form: they are invariant to rotation and translation, but contain information on shape and overall scale. Another set of form variables would be the distance between A and B, along with the two perpendicular distances d and e. Angles are often avoided in this context because of their demanding statistical properties.

Fig. 1
figure 1

Three landmarks A, B, and C constituting a triangle

Shape variables, in contrast, are independent of the overall size of a configuration and are often constructed by dividing some interlandmark distances by a measure of scale. In the absence of shape differences, overall size or scale is an unambiguous concept and may be measured, e.g., by any arbitrary interlandmark distance. When shape varies across the observed geometries, size is less clearly defined and requires either the (arbitrary) choice of a representative single distance measure (usually one that is less affected by the observed shape variation) or a quantity based on a series of single size measures (such as Centroid Size, see below, or the average of several distance measurements). Another method to describe the shape of a triangle is called two-point shape coordinates or Bookstein shape coordinates (Bookstein 1991). Translate, scale, and rotate the triangle in Fig. 1 until A has coordinates (0,0) and B (1,0). The ensuing two coordinates of C (the shape coordinates) sufficiently describe the shape of the triangle. Alternatively, they can be computed as d and e divided by the distance AB (the baseline size). A more general way to compute shape coordinates, the Procrustes superimposition, is described below.

The shape of a triangle can be captured using two variables only and thus has two degrees of freedom (df) whereas the description of its form requires three df. In addition to the information about the triangle’s shape, the six coordinates of the three landmarks contain information on scale (1 df), orientation (1 df), and position (2 df)—so called nuisance parameters. In general, for p landmarks in k dimensions, shape has pk-k-k(k-1)/2-1 df, that is four (for two-dimensional data) or seven (3D) df less than the number of landmark coordinates.

The two dimensions of shape for a triangle constitute a two-dimensional shape space, in which a single point represents the shape of a triangle. Interestingly, this shape space is not a flat plane but has the form of the surface of a sphere (Fig. 2); it is called Kendall shape space after the Scottish mathematician who discovered it (Kendall 1981, 1984; Slice 2001). For more than three landmarks Kendall shape space is a more complex Riemannian manifold (a generalization to higher dimensions of a curved surface in three dimensions). Kendall (1984) demonstrated that if the vertices of a shape are independently and identically distributed spherical normal variables, then the distribution of shape is isotropic (the same in every direction) in Kendall’s shape space (see also Bookstein 1991; Dryden and Mardia 1998). This is an important property, guaranteeing that isotropic measurement error does not induce an apparent “signal” in shape space.

Fig. 2
figure 2

a The vertices of one hundred random triangles (independently and identically distributed normal variables) with an equilateral triangle as the mean form. b The shapes of these triangles are isotropically distributed in Kendall’s shape space, which has the form of the surface of a (hemi)sphere

An empirical analysis based on linear statistical methods may thus be confounded by the non-linearity of shape space, i.e., by variation in the pk shape coordinates due to the curvature of shape space rather than due to actual shape differences among the observed forms. Because variation in biological shape is relatively small even when observed across a wide range of different taxa, it is possible to make a good linear approximation to the non-linear shape space (Rohlf 1999; Marcus et al. 2000). This linear space is of the same dimensionality as shape space and can be viewed as tangent to it, where the point of tangency is at the reference shape (usually taken as the sample average). Euclidean distances in this tangent space closely resemble the distances in Kendall’s shape space, so that the shapes projected into tangent space can be used for analysis with standard multivariate methods. Furthermore, the Euclidean geometry of tangent space gives rise to meaningful notions of distance, length, and direction, which are deeply entrenched in the rhetoric of modern evolutionary and theoretical biology. For example, distances (i.e., similarities) among individual shapes or group mean shapes can be mutually related, and the length and direction of developmental or evolutionary trajectories can be meaningfully compared. This is not the case for many classical “morphospaces” (Mitteroecker and Huttegger, in press).

There are two principal ways to construct a tangent space, but because of its geometric properties an orthogonal projection onto the space perpendicular to the vector of shape coordinates of the reference shape is preferred in most cases (see Rohlf 1999, for more details). Many practitioners of morphometrics simply ignore the curvature of shape space and apply standard multivariate techniques to shape variables without prior projection into tangent space. This may not be problematic for many data sets (most linear methods, such as principal component analysis, linearly approximate the curved shape space) but may lead to serious artifacts when the variation within the data is largely constrained to a single direction of shape space (e.g., a single growth trajectory; see Mitteroecker 2007).

Procrustes Superimposition

For most practical applications, the parameters describing the shapes for a sample of homologous landmark configurations are estimated by a Procrustes superimposition. This procedure is a least-squares oriented approach involving three steps (Fig. 3):

Fig. 3
figure 3

The three steps of Procrustes superimposition: translation to a common origin, scaling to unit centroid size, and rotation to minimize the sum of squared Euclidean distances among the homologous landmarks. The resulting landmark coordinates are called Procrustes shape coordinates

  1. 1.

    Translation of the landmark configurations of all objects so that they share the same centroid (the coordinate-wise average of the landmarks of one form). Usually, this common centroid is sent to the origin of the coordinate system.

  2. 2.

    Scaling of the landmark configurations so that they all have the same Centroid Size—the square root of the summed squared deviations of the coordinates from their centroid. Centroid Size is a measure of scale for landmark configurations, which has been shown to be approximately uncorrelated with shape for small isotropic landmark variation (Bookstein 1991; Dryden and Mardia 1998). As a convention, Centroid Size is set to one for all landmark configurations.

  3. 3.

    When superimposing two landmark configurations, one of the two centered and scaled configurations is rotated around its centroid until the sum of the squared Euclidian distances between the homologous landmarks is minimal. For more than two forms, this algorithm has been extended to Generalized Procrustes Analysis (GPA), where the rotation step is an iterative algorithm (Gower 1975; Rohlf and Slice 1990). First, all centered and scaled landmark configurations are rotated to one arbitrary configuration of the sample using the same least squares approach as above. Subsequently, the resulting coordinates are averaged and all configurations are rotated again to fit this tentative consensus. These new coordinates are averaged as an updated template for the next iteration. The algorithm is repeated until convergence, which is usually reached after a few iterations.

The coordinates of the resulting centered, scaled, and rotated landmarks are called Procrustes shape coordinates; their average (the consensus configuration) is the shape whose sum of squared distances to the other shapes is minimal and is thus the maximum likelihood estimate of the mean for certain statistical models (Dryden and Mardia 1998). The individual differences from the average shape are called Procrustes residuals. The Euclidian distance between two sets of Procrustes shape coordinates (i.e., the square root of summed squared coordinate-wise differences) is referred to as Procrustes distance and denotes the (dis)similarity in shape between two landmark configurations.

Even though scaling the specimens to unit Centroid Size resembles several other standard approaches to size correction, it is not the actual least squares solution. The algorithm described above is therefore sometimes referred to as partial Procrustes fitting. A smaller sum of squared deviations among the landmark configurations can be achieved by constraining the size of a configuration to the cosine of the angle between the vector of shape coordinates of that specimen and the vector of the mean shape. This scaling approach results in the so-called full Procrustes fit (Rohlf 1999). In most published analyses, however, only partial Procrustes fitting is considered. Recently, Theobald and Wuttke (2006) proposed a maximum-likelihood approach to Procrustes superimposition (but see also Bookstein 2007).

Due to the scaling step in the course of Procrustes superimposition, the shape coordinates do not possess any information on overall size—statistical results can only be interpreted in terms of relative sizes within the observed specimens. Yet, in many biological studies overall size is an integral part of the assessed morphology and the analysis should be carried out in form space rather than in shape space. In Mitteroecker et al. (2004a, 2005) we have demonstrated that the Procrustes shape coordinates can be augmented by the natural logarithm of Centroid Size to construct Procrustes form space (also called size-shape space). As for shape space, independently and identically distributed spherical normal variation (such as measurement error) around a mean form results in an isotropic distribution in form space. Form space should be used when both size and shape is of scientific concern, like in classification studies or the analysis of development and allometry. Developmental or evolutionary dissociation of size and shape, such as in the context of heterochrony (Zelditch 2001), can best be studied when contrasting results in shape space and form space (Mitteroecker et al. 2004b, 2005; Gerber et al. 2007).

Deformation Grids

Visualization of shape differences and shape changes is a primary aim and major strength of geometric morphometrics. Morphological comparisons by deformations grids date back to at least Albrecht Dürer’s 1528 book “Vier Bücher von menschlicher Proportion” and D’Arcy Thompson’s (1915, 1917) pictorial approach of “Cartesian transformations”. Based on a mapping of homologous point locations between forms, Thompson constructed deformation grids illustrating how a part of one organism may be described as a distortion of the same part in another individual. Dürer’s and Thompson’s approaches are visually appealing but their drawings were still made by hand without any reference to a formal algorithm. It required fast personal computers and the appropriate mathematics to apply the idea of deformation grids on a computational basis. After several other, meanwhile discarded attempts, Bookstein (1989, 1991) introduced the method of thin-plate spline (TPS) interpolation to compute deformation grids in the style of Thompson and Dürer. This algorithm, borrowed from material physics, computes a mapping function between two point configurations that maps the measured points exactly while the space in-between is smoothly interpolated. The notion of smoothness is approached by minimizing the bending energy of the deformation, a scalar quantity computed as the integral of the squared second derivatives of that deformation. The TPS formalism is also central to the semilandmark algorithm and the estimation of missing data in morphometrics (see below).

The TPS interpolation function from a template configuration to a target configuration is usually applied to the vertices of a regular grid so that the shape differences between the two geometries can be read from the deformation of this grid (Figs. 4 and 5a). When the actual shape differences are subtle, the deformation can be extrapolated by an arbitrary factor to ease the interpretation of the grid. In the course of computation, the thin-plate spline function is applied to each coordinate axis separately and can thus be used for both two-dimensional and three-dimensional data. Deformation grids have proved less effective for visualizing 3D shape differences, but because the mapping function can be applied to any points in the vicinity of the template landmarks, the algorithm can also be used to deform a surface model (the vertices of a surface triangulation) of the template specimen. A sequence of warped surfaces can provide a useful alternative to deformation grids for describing three-dimensional shape and form differences (Fig. 5b).

Fig. 4
figure 4

A template configuration (left) and a target configuration (right) of five landmarks each. The deformation grid on the right illustrates the thin-plate spline function between these configurations as applied to the left regular grid—it is a visualization of the differences between the two shapes

Fig. 5
figure 5

a In a sample of 19 landmarks digitized on each of 691 cichlid fish, the first principal component is visualized by a thin-plate spline (TPS) function, applied to both a deformation grid and an image of a fish. The left image represents a warp from the consensus configuration (the average) to the negative direction of the principal component whereas the right image is the opposite deformation, thus representing a warp towards the positive direction (from Herler et al., in press). b Visualization of the average postnatal growth in the hominoid cranium based on a shape regression on log Centroid size in a sample of 96 landmarks and semilandmarks measured on 268 specimens. Starting from an average newborn, the TPS function is applied at four equal steps along the vector of average growth to the vertices of a surface triangulation, which was extracted from a CT-scan (from Mitteroecker et al., 2004a). Note that while the TPS function is applied directly to the deformation grid and the surface triangulation (forward warping), the fish image is deformed by an unwarping (see main text)

The TPS function can also be applied to the pixels of an image or to the voxels of volumetric data as derived from CT or MRI scans. However, when warping the pixel locations from the template to the target space according to the two landmark configurations, pixels may overlap in the target image (in areas of compression) or positions in the image may be empty (in areas of expansion). To avoid such fragmented images, the TPS algorithm is often used instead to unwarp an image. The pixel positions of the target image are warped to the template, identifying the pixels that correspond across the two images. The gray values or the color values of the target pixels are then substituted by the corresponding values in the template image. The resulting unwarped image is close (but not identical) to a respective forward warping and is a continuous image with no gaps or wholes (Fig. 5a).

In principle, bending energy is a measure of shape difference between two landmark configurations, which does not require Procrustes superimposition of the configurations. It gives strong weight to shape differences at a small geometric scale, while not taking into account affine shape differences at all (shape deformations with infinite scale). But bending energy is not a metric measure of distance (e.g., it is not symmetric) and hence is usually not used for statistical analysis. The TPS formalism can also be applied to decompose shape deformations into a range of geometrically independent components (partial warps) with different geometric scale and hence different bending energy (Bookstein 1989, 1991). A decomposition of the mean form (principal warps) can be used as an orthogonal basis to span tangent space, but is otherwise of limited biological relevance (Rohlf 1998; Monteiro 2000).

Semilandmarks

All geometric morphometric tools strictly require that the digitized landmarks are homologous across specimens, i.e., they represent the same biological locations in every individual. Thus, anatomical landmarks need to be well defined in all two or three dimensions of the digitizing space and must be identifiable on all specimens in the sample. Typical landmark locations are small cusps, invaginations, or intersections of tissues such as bony sutures. In many applications, however, point locations that fulfill all these requirements are not evenly distributed across a biological organism. A vertebrate skull, for example, has a large number of anatomical landmarks in the face but only a few points can be defined unambiguously on the smooth braincase. Due to this limitation, most geometric morphometric analyses of skull form have focused on the face and the cranial base. The concept of semilandmarks (also called sliding landmarks), first introduced in an appendix to the “Orange Book” (Bookstein 1991), was invented to extend landmark-based statistics to smooth curves and surfaces. It was first applied to two-dimensional outlines (Bookstein 1997) in a study of shape variability of sections of the human corpus callosum (the structure that connects the left and right hemispheres of our brain). In Gunz (2001, 2005) and Gunz et al. (2005) we extended the algebra to curves and surfaces in three dimensions. Alternative computational approaches to semilandmarks can be found in Frost et al. (2003) and Reddy et al. (2005).

Even though it is impossible to define an exact point-to-point correspondence between landmarks on smooth curves or surfaces, the semilandmark algorithm requires that the structures themselves are homologous. While the location of a landmark along the curvature may not be identifiable, its coordinate perpendicular to the curvature is well defined and driven by the observed morphology. The semilandmark algorithm discards information derived from the arbitrary spacing of semilandmarks along the curves or surfaces, thus identifying the coordinates of biological interest—the spacing of semilandmarks is produced as a byproduct of the statistical analysis itself. This is achieved by allowing the points to slide along their curve or surface until some measure of shape difference among the configurations is minimized. The two main semilandmark algorithms differ in their optimization criteria: in one case the bending energy of the thin-plate spline between each specimen and the sample mean shape is minimized whereas the Procrustes distances are minimized in the other version of the algorithm.

To linearize the minimization problem, the semilandmarks do not slide on the actual curve or surface but along the tangent vectors to the curve or the tangent planes to the surface (Fig. 6). Sliding of semilandmarks is an iterative process repeating the following three steps: (1) computation of a tentative sample mean shape and of the tangent vectors for each semilandmark in each specimen, (2) sliding of semilandmarks along the tangents to minimize either the Procrustes distances or the bending energy to the mean, and (3) placing back the slid semilandmarks to the nearest point on the curve or the surface, as they may slip off the curvature (for relatively smooth curves this third step may be omitted).

Fig. 6
figure 6

When the location of a landmark on a smooth curve or surface cannot be identified clearly, it may be treated as a semilandmark that is allowed to slide along its curvature—only the position perpendicular to the curve or surface carries a biological signal. As a linearization of the underlying algorithm, the landmarks do not slide along the actual curvature but on the tangent structures. For curves, such as a section through a human corpus callosum in a or three-dimensional ridge curves on a human cranium in b, the semilandmarks are constrained to tangent vectors (the black and gray straight lines, respectively). Semilandmarks on surfaces, such as the cranial vault in b, are allowed to slide on tangent planes defined by two vectors per semilandmark

Semilandmarks have to be equal in number across the sample and their starting positions should be in rough geometrical correspondence. Clearly observable curves on surfaces, such as ridges, should be treated as curves instead of surface points. In general, the sampling of semilandmarks depends on the complexity of curves or surfaces and on the spatial scale of shape variation that is of interest. Sampling experiments can help finding an “optimal” number of semilandmarks in the sense of how much information additional points would contribute (Katina et al. 2007).

To arrive at the same number of semilandmarks in the same order on each specimen, it is convenient to begin with points equidistantly spaced along outline arcs, e.g., through automatic resampling of a polygonal approximation to the curve. Techniques for surfaces differ substantially from those for curves in that, except for planes and cylinders, there is no straightforward analogue to the notion of “equal spacing”. We recommend beginning with a hugely redundant sample of points on each surface, e.g., point clouds generated by a surface scanner or extracted surfaces from CT or MRI data. Alternatively, they can be produced just by “scribbling” around the surface with a digitizing device such as a Microscribe. On a single reference specimen, one then produces a mesh of far fewer and relatively evenly spaced points by thinning the redundant point cloud. Points on this mesh should be denser near ridges of the surface, especially if they are not treated as curves. This reference mesh is then warped into the vicinity of the surface of every other specimen using a thin-plate spline interpolation function based on the anatomical landmarks and the (yet unslid) curve-semilandmarks. On the surface representation of each specimen the points nearest to the warped mesh are taken as starting positions for the semilandmark algorithm. Through sliding the semilandmarks acquire geometric homology or correspondence and can be used in the subsequent statistical analyses just as if they were classical landmarks.

Statistical Analysis

Statistical analysis in geometric morphometrics is performed using the Procrustes shape coordinates (or other representations of the organism’s geometry) and differs in several respects from traditional approaches. In traditional morphometrics, multivariate statistics is applied to a set of measurements including distances, distance ratios, angles, volumes, counts etc. As such variables often do not share common units and comparable ranges of variation, they typically are log transformed and standardized to unit variance. Statistical analysis is thus usually based on correlation matrices and comprises the full range of multivariate techniques and statistical tests (see, e.g., Blackith and Reyment 1971; Marcus 1990; Oxnard 1983). In geometric morphometrics the shape variables all possess the same units so that analyses are based on the covariance matrix and there exists a well-defined metric (the Procrustes metric). Results of multivariate methods that preserve this metric, such as principal component analysis, multivariate regression, and partial least squares, can be visualized as actual shapes or shape deformations in the geometry of the original specimens. The algebra of the statistical methods is the same for two- and three-dimensional shape data; only the kind of visualization may differ (see above).

Some methods which are popular in traditional morphometrics, such as multiple regression, canonical correlation analysis, and canonical variate analysis are not equally useful in geometric morphometrics for several reasons: most importantly they do not preserve the Procrustes metric, discarding one of the core strengths of geometric morphometrics. In addition, their computation may be unreliable or impossible because these methods all involve the inversion of the landmark variance-covariance matrix. To compute a reliable matrix inverse, the number of cases must exceed the number of variables, which is often not the case in modern morphometrics. Also, Procrustes shape coordinates usually are highly correlated (especially when using semilandmarks), which results in a singular covariance matrix that cannot be inverted. Even when the number of variables is reduced, such as by principal component analysis (see below), the results of these methods may be difficult to interpret and to visualize directly as shape deformations.

For example, canonical variate analysis (CVA), a popular method originally designed to classify unknown individuals into pre-specified groups, is frequently used in morphometrics in the style of an ordination technique and to assess whether two or more groups differ in their mean tendency. But when the number of variables is close to the number of cases—a common situation in geometric morphometrics—CVA will always separate groups even if they actually possess the same mean. As a consequence, such approaches are best avoided in geometric morphometrics and a smaller range of multivariate methods has been established as core statistical techniques.

Principal component analysis (PCA) is a method to reduce a large set of variables to a few dimensions that represent most of the variation in the data. PCA is computed by an eigendecomposition of the sample covariance matrix and is a rigid rotation of the data preserving the Procrustes distances among the specimens. Principal component scores are the projections of the shapes onto the low-dimensional space spanned by the eigenvectors. They can be plotted as two- or three-dimensional graphics and allow one to assess group differences, growth trends, outliers, etc., in the data without incorporating prior information such as group affiliation; they are based on the shape or form variables only. The eigenvectors or principal components contain the weightings for the linear combinations of the original variables and can be visualized as actual shape deformations (relative warps, Bookstein 1991; see also Fig. 5a). However, the principal components are statistical artefacts largely depending on the composition of the actual sample, and should not be interpreted as representing biologically meaningful factors (cf. Mitteroecker et al. 2005; Mitteroecker and Bookstein 2007; Adams et al., in review). There is one exception to this rule in morphometrics: the first principal component of a single species or population often represents allometry, the shape variation induced by overall size variation (Klingenberg 1998). But this is only the case if allometric variation is the dominant factor in the data, such as in an ontogenetic sample, otherwise multivariate regression should be preferred to estimate allometry.

Multivariate regression of the shape variables on an external variable (shape regression) is the method of choice to assess the influence of a single factor such as a functional or environmental variable, age, or the dose of a drug on shape (Bookstein 1991; Monteiro 1999). The regression of shape on the logarithm of Centroid Size is the optimal measure of allometry (Mitteroecker et al. 2004a). Multivariate regression is not sensitive to the number of dependent shape variables or to their covariance structure and the resulting vector of regression coefficients (quantifying the average effect on shape) can be visualized as shape deformation (Fig. 5b).

Partial least squares (PLS) is a method to assess relationships among two or more blocks of variables measured on the same entities (Wold 1966; Rohlf and Corti 2000; Bookstein et al. 2003). PLS yields linear combinations that optimally (in a least-squares sense) describe the covariances among the sets of variables and so provide a low-dimensional basis to compare different blocks of variables. Unlike canonical correlation analysis, the computation does not involve the inversion of a matrix and the results can be interpreted conveniently in terms of path models and can be visualized as shape deformations (singular warps). When all blocks of variables are comprised of shape coordinates, PLS may be used to assess the “morphological integration” among the respective anatomical components (Bookstein et al. 2003; Bastir and Rosas 2006; Gunz and Harvati 2007). As PLS only gives the patterns of shape variation with the highest mutual covariance, but not the corresponding amounts of these deformations (the singular vectors are all scaled to unit length), a separate scaling step may be necessary (Mitteroecker and Bookstein 2007, 2008). But the blocks may also be comprised of other variables, such as functional, environmental, or behavioral measures (e.g., Bookstein et al. 2002; Frost et al. 2003, Manfreda et al. 2006), and PLS can be used to identify the latent variables underlying the association among shape and those factors.

Many geometric morphometric analyses are based on randomization tests, such as permutation or bootstrap tests (Manly 1991; Good 2000), rather than parametric methods to assess the statistical significance level of a given hypothesis. Most parametric tests require more cases than variables and a specific (usually a normal) distribution of the variables. Randomization tests, in contrast, are free from these restrictions as long as the cases are sampled independently. Furthermore, test statistics can be designed even for complex hypotheses and compared with their permutation distribution.

In summary, statistical analysis in geometric morphometrics preserves the original geometry of the landmark configurations and hence permits the visualization of results as shape deformations. Together with the relatively large amount of shape variables, this gives rise to a specific exploratory style of analysis. Previously unknown shape features can be inferred from visualizations such as TPS deformation grids (e.g., Bookstein 2000) and vague notions of a signal can be “calibrated” (Martens and Naes 1989) with the present data. Instead of interpreting differences along single shape coordinates or principal components, it is more powerful to assess actual group mean differences, shape regressions, and singular warps, which can be visualized and directly related to the experimental design. Traditional morphometrics, in contrast, is based on some selected measures and the results are naturally restricted to those few variables. Such classic multivariate analyses may be sufficiently effective in assessing shape features known prior to the data collection but are unlikely to produce novel and unexpected findings.

Asymmetry

For bilateral objects, shape space can be partitioned into a subspace of symmetric shape variation and a subspace of asymmetric variation. The actual computation is based on a Procrustes superimposition of the landmark configurations together with their relabeled reflections (Bookstein 1991; Klingenberg and McIntyre 1998; Mardia et al. 2000). The Procrustes distance between a shape and its reflection is a measure of asymmetry; it is zero only for perfectly symmetric shapes. A shape may be “symmetrized” by averaging the shape and its superimposed reflection (Fig. 7). Constraining an analysis to the symmetric part of shape variation reduces the degrees of freedom by discarding asymmetric shape variation, including asymmetric measurement error, which may not be of interest for the analysis.

Fig. 7
figure 7

Analysis of asymmetry in geometric morphometrics: a an asymmetric configuration of 64 landmarks digitized on a photograph of a human face. The configuration in b is the average of the asymmetric shape a and its reflection c after Procrustes fitting, resulting in a perfectly symmetric shape. The Procrustes distance between a and c, which equals two times the distance between a and b, is a measure of total asymmetry. The deformation grid d from the symmetric consensus to the asymmetric shape in a visualizes the pattern of asymmetry. This deformation is extrapolated by a factor of three to ease its interpretation

In a sample of shapes, asymmetric variation can be decomposed further into directional asymmetry (the asymmetry of the population mean) and fluctuating asymmetry (asymmetric deviations from the population mean), based on a multivariate analysis of variance. Fluctuating asymmetry has been used to assess developmental instability (“intra-individual variation”) in many recent investigations (e.g., Debat et al. 2000; Hallgrímsson et al. 2003; Schaefer et al. 2006; Willmore et al. 2006).

Missing Data in Morphometrics

From a statistical point of view, landmarks that cannot be measured because anatomical structures are broken off or deformed are termed “missing data”. As geometric morphometric methods require all measured specimens to have the same number of landmarks at homologous positions, the variables or the cases with missing values have to be either excluded from the analysis, or the missing values need to be estimated based on the available data. There exists a large body of literature on methods for estimating population parameters such as means or regressions in the presence of missing values (for reviews and examples see Schafer 1997; McLachlan and Krishnan 1997; Schafer and Graham 2002). In morphometrics, particularly when applied to the paleo-sciences, the estimates of the missing values may be of interest themselves, e.g., as a way to reconstruct fragmented fossils. In Gunz (2005) and Gunz et al. (2004, in review) we have laid out a framework for estimating missing shape coordinates, exploiting the high redundancy of shape variables—especially when the measurement points are closely spaced—and taking into account prior biological knowledge about functional and geometric constraints, symmetry, and morphological integration.

The easiest yet most powerful way to deal with missing and deformed anatomical structures is to restore bilateral symmetry. When parts are missing on both sides or in the symmetry plane, one can predict the coordinates of the missing landmarks using information from complete cases. Statistical reconstruction is aimed at optimizing the statistical likelihood of the estimate, using an extension of a common missing data algorithm, the expectation-maximization method (Dempster et al. 1977), for Procrustes shape variables. Geometric reconstruction employs the smoothness properties of the thin-plate spline to estimate the missing coordinates based on a single reference configuration. A TPS interpolation is computed based on the subset of landmarks that are available in both the reference and the incomplete specimen. This interpolation function is used to map the missing landmarks from the reference onto the incomplete target, placing the landmark estimates so that the deformation between the reference and the incomplete specimen is as smooth as possible (minimum bending energy). The reference configuration can be an actual specimen or a Procrustes group average and may resemble the incomplete specimens in properties such as group affiliation, age, or geographic origin. Naturally, choosing different references will result in slightly different reconstructions. Whether similar conclusions may follow from a variety of realistic alternative reconstructions is testable in the course of the analysis: the distribution of differently reconstructed shapes reflects the uncertainty due to the missing data and the sensitivity to prior assumptions.

Morphometric Resources in the Internet

There emerged a range of free software for geometric morphometric analysis in the last two decades. The Stony Brook morphometrics site (life.bio.sunysb.edu/morph) provides links for most of these software packages along with morphometric data sets and bibliographies. The TPS program series by James Rohlf, offering possibilities for digitizing landmark, Procrustes superimposition including semilandmarks, deformation grids, image warping, and shape regression, became the standard tool for two-dimensional landmark data. Several programs have been developed to handle three-dimensional landmark data, which mainly requires more extensive visualization tools. The most comprehensive free software packages are Morpheus et al. by Dennis Slice, Morphologika by Paul O’Higgins, IMP by David Sheets, Past by Øyvind Hammer, MorphoJ by Chris Klingenberg, and APS by Xavier Penin. Morpheus et al., which offers elaborate 3D visualizations, and MorphoJ are written in Java and thus are platform-independent. The Edgewarp3D package from Bill Green and Fred Bookstein runs under Linux and Mac OS X and provides tools for digitizing landmarks on 3D surfaces and volumes, including semilandmarks on curves and surfaces, together with many visualization options. With the software package Landmark Editor (developed by UC Davis in collaboration with the NYCEP morphometrics group) one can measure landmarks on 3D surfaces and capture semilandmarks by creating “patches” with equal point count. Ian Dryden programmed several procedures for statistical shape analysis in R (library “shapes”) and David Polly as well as Philipp Gunz and Philipp Mitteroecker developed a series of functions for geometric morphometrics in Mathematica.

Further useful links can be found at www.morphometrics.org by Dennis Slice, including a possibility to subscribe to the morphmet mailing list. An extensive morphometrics glossary by Slice and co-workers is available at life.bio.sunysb.edu/morph/glossary/gloss1.html.