19. Shape Models II: The Thin Plate Spline
Written by Norm MacLeod - The Natural History Museum, London, UK (email: n.macleod@nhm.ac.uk). This article first appeared in the Nº 73 edition of Palaeontology Newletter.
Introduction
In the last column we developed a simple method for expressing the results of complex, multivariate, geometric morphometric ordinations as form and/or shape models. These were configurations of landmark points that exist at discrete coordinate locations in the multivariate linear spaces we typically use in geometric morphometrics to portray similarities and differences in the form and/or shape of a sample of specimens. In this column we’re going to develop a more mathematically sophisticated way to do the same thing that, at least on a superficial level, makes connection with the ‘deformation grid’ approach to shape modelling developed by D’Arcy Thompson in his classic treatise On Growth and Form (1917, 1942). But before we get into the mathematics of this, now conventional, approach to expressing the results of a comparison between forms/shapes, and to avoid later confusion regarding the degree to which this convention realizes the Thompsonian ideal, it’s worth taking a moment to review what Thompsonian transformation grids are and why they were developed.
In his original (1917), and in the later, expanded (1942) edition, Thompson’s goal was to “correlate with mathematical statement and physical law certain of the simpler outward phenomena of organic growth and structure or form while all the while regarding the fabric of the organism, ex hypothesi, as a material and mechanical configuration” (Thompson 1917, p. 17). In other words, Thompson sought to “see how, in some cases at least, the forms of living things, and of parts of living things, can be explained by physical considerations and to realize that in general no organic forms exist save such as are in conformity with physical and mathematical laws” (ibid, p. 15). Thompson was, of course, aware of evolutionary theory and agreed that natural selection operated to sort mechanically efficient from inefficient forms in the manner Darwin had suggested. But he bridled at the idea that every aspect of a form is now, and always has been, under direct adaptive scrutiny, preferring to believe that some aspects of form owe their origin to the physical forces with which they must contend.1 Thompson saw organic form as a ‘diagram of forces’ from which inferences can be made regarding the nature of the forces that act upon it now or that have acted upon it in the past. Using this force metaphor, Thompson saw the mathematical comparison of forms as a way of deducing how these fields of forces changed during both ontogenetic and evolutionary history.
Thompson’s proposed method of force-field analysis was to take two simple line drawings of species’ bodies or some corresponding structural element therefrom (e.g., copepod, ungulate cannon bone, leaf). For convenience we’ll label one form as the ‘reference’ and the other as the ‘target’. In order to better visualize the nature of the geometric transformation Thompson superimposed a rectilinear grid on the reference form. He then worked out simple sets of mathematical transformations that would map point location coordinates of the reference into topologically corresponding locations on the target. Applying those same mathematical transformations to the coordinate locations of the grid vertices Thompson obtained a striking image-based summary of the implied geometric transformation (Fig. 1). In this view Thompson perhaps reflects the same level of discomfort with hyperadaptationist arguments criticized more recently by Gould and Lewontin (1979), among others.
Figure 1. Example Thompsonian transformation grids specifying uniform (upper) and non-uniform (lower) transformation functions. For each comparison the reference form is located on the left and the target form on the right. Redrawn from Thompson 1917.
It is clear from Thompson’s many statements throughout the last chapter of On Growth and Form that he regarded the mathematical transformation as pertaining to, and being constrained by, all mathematical points comprising the line drawing and that he respected the principle that biological homology pertained to structures, but not necessarily individual point locations on structures. Rather, it was the configuration of the entire ensemble of mathematical points — represented diagrammatically by the superimposed grid — that he looked to in judging whether he had devised biologically reasonable formulae for a particular form transformation. Similarly, it is clear the only purpose served by the mathematical grid was to passively express the overall geometry of the transformation in the manner of a deformed, map-like coordinate system.
As illustrated in Figure 1, Thompson used his approach to provide examples of both linear (uniform) and non-linear (non-uniform) transformation modes. Like the later ‘relative growth’ studies of Otto Snell, Julian Huxley and Georges Teissier (see below), the thing that impressed morphologists about Thompson’s transformation grids was the fact that seeming complex form changes appeared to be able to be described accurately by simple mathematical transformations applied consistently to all point locations over a form. This suggested to many at the time that the underlying principles and/or determinants of morphological change might be simple when expressed in, or studied using, the language of mathematics.
While Thompson’s transformation-grid approach resulted in the creation of compelling diagrams — so much so both his original set of drawings and many subsequent variations of them have been reproduced in countless books on biology and evolution despite the fact that the physicalforce theory these drawings represent is almost never discussed in those same texts — his geometric approach to the analysis of form never caught on during his lifetime. Thompson himself provided some guidance regarding how to operationalize his transformation grids, which he thought of as a visual tool akin to a modern-day spatial morphing algorithm. Those algorithms subdivide an image into a set of points and then smoothly map a subset of these between a reference and target form with their difference displacements informing the displacement of intermediate points via simple linear interpolation. For example, it is this linear interpolation approach to transformation grid analysis that Thompson used to create his morphlike model of the complex geometric transition between Hyracotherium and modern-day Equus (Fig. 2).
Figure 2. Thompson’s transformation grid analysis for the transition from the Eocene Hyracotherium skull form (A) to the modern Equus skull form (B). Note the representation of hypothetical intermediate stages of the transformation via linear interpolation (C). From Thompson 1917.
Although Thompson’s original, interpolation-based approach to the realization of transformation grids was fine if all you wanted to do was map one form into another, it was not well suited to the summarization of geometric information across a larger sample of data. Arguably Huxley (1924) and Teissier (1929, see also Snell 1892 and Huxley 1932) were more successful in developing an analytic approach to the general problem of form variation than was Thompson. But the logistic regression equations used by students of relative growth — or allometry in modern parlance — were not used to create graphic models of form change with anything like the visual impact of Thompson’s grids.
Curiously, this failure to capitalize on the modelling capabilities of regression-based methods when treating morphological data remained in place for a half-century during which time Huxley and Teissier’s bivariate regression-based approach to form analysis was extended to the multivariate case (via PCA, see Jolicouer and Mosimann 1960) and holistic modelling approaches were developed for other aspects of morphological analysis (see Olson and Miller 1958). As we have seen in the last column, all the mathematical machinery for implementing at least some aspects of useful geometric shape modelling was in place by the 1960s. Yet, no new developments in this area took hold until 1980s despite a few attempts to formulate an explicitly Thompsonian modelling approach in the form of morphological trend surfaces (Sneath 1967) and biorthogonal grids (Bookstein 1978).
In retrospect there appear to be two reasons for this. The first was that the largest school of morphometrics (multivariate morphometrics) tended to look to the communities of statisticians and psychometricians for methodological guidance, neither of which were particularly interested in creating morphological models. The second was that, ever since the 1930s, the tradition in bivariate and multivariate morphological studies was to analyze pairs or sets of distances between landmark locations rather than configurations of Cartesian coordinate locations scattered over a sample of forms. Once the power of shape coordinates had been established by Bookstein (1986) and the outlines of shape theory had begun to emerge (see Kendall 1984), the stage was set to renew the search for an analytic method that could combine the intuitive appeal of Thompson’s transformation grids with the equally popular, and far more powerful, tools of multivariate morphometrics. The key insight that allowed this new approach to shape modelling to be realized was specification of a new spatial metaphor for shape similarity.
In traditional multivariate analysis the similarity between two objects is quantified by calculating the distance between them across all variables. This is fine for a quick-and-dirty summary of form differences, but lacks the ability to track varying patterns of size/shape similarity and difference in different regions of the forms. In the days when morphometricians characterized forms using sets of linear distances between landmarks, this distance-based metaphor seemed both natural and practical. After all, distances are simply magnitudes. There is no information about geometry in a list of distance values. Since geometry can’t be reconstructed precisely from a table of distance values there wasn’t any point in worrying about shape models. But with the move to characterizing forms using the coordinate values of the landmarks themselves — and especially the transformation of landmark coordinate values to Procrustes shape coordinate values via standardization for position, scale, and rotation — it became possible to represent the form and shape similarities or differences between any two objects precisely in a manner that retained the fundamental geometry of the landmark configurations.
Figure 3. Stages in the landmark-based comparison between shapes. (A) Cranidial landmark configuration for the trilobite genus Acaste. (B) Cranidial landmark configuration for the trilobite genus Calymene. (C) Procrustes superposition of the Acaste (red) and Calymene (green) shape coordinates. (D) Shape difference vectors between corresponding reference (red) and target (green) shape coordinate locations. See text for discussion.
Imagine two forms defined by corresponding sets of landmark coordinates. As above, we’ll call one the reference form and the other the target (Fig. 3A-B). Procrustes superposition of these landmark sets transforms the forms into shapes and brings corresponding landmarks into positions of maximal correspondence (= minimal sum of squared deviations, Fig. 3C). The differences between the reference and target shapes can be visualized as a set of difference vectors between the reference and target shapes at each coordinate location (Fig. 3D). Now, rather than summing these differences up to produce a statistical estimate of shape-difference as we’d do in traditional multivariate morphometrics, let’s select the reference shape and express the difference between it and the target shape as a set of vectors with the same displacement as in Fig. 3D, but rotated such that the difference vectors are parallel to the z-axis of a three-dimensional coordinate system. In this system the z-axis expresses the shape difference between the reference and target. The resulting figure (Fig. 4A) expresses the shape difference as a set of stalks at each landmark location. Either the reference or the target shape can be used as the basal shape for constructing this sort of shape-difference diagram. Note that the geometry of the landmark configuration is retained in this graphic expression of shape difference and that we can easily identify which regions of the two shapes are more similar to, or more different from, one another.
Figure 4. Shape-difference diagrams between the Acaste (reference) and Calymene (target) cranidial landmark configurations. (A) Rotation of the reference-to-target shape-difference vectors shown in Fig. 3D to the z-axis to form a three-dimensional representation of shape difference. (B) Spline-estimated surface of the shape-displacement data with colors representing regions of differential curvature. (C) Same spline surface as in B, but this time with the surface represented as a deformed rectilinear grid. In all three diagrams the reference configuration was used as the basis for the shape-difference graphic.
The surfaces shown in Figures 4B and 4C are standard parametric cubic splines. These aren’t terribly useful for summarizing shape difference because they specify elastic deformations in which the shape of the underlying mesh, when viewed from the direction of the shape-displacement (z) axis, is held constant. This amounts to an isomorphic projection of the 3D mesh onto the x,y plane. Instead, the surface interpolation method of choice among morphometricians is an advanced type of polyharmonic spline called a thin plate spline (TPS, see Duchon 1976, 1977). The TPS attempts to mimic the behavior of a defect-free, uniform, and infinitely thin metal plate that is bent in the z-direction to conform to the geometry of the shapedisplacement vectors. This metal-plate metaphor is important because, unlike elastic surfaces, metal plates bend in ways that minimize the energy required to achieve the bend in all directions over the entire plate. In applying this physical metaphor to shape analysis, the TPS represents the surface of minimal bending energy implied by the transformation of one shape into another.
Before we get into the equations for calculating a TPS surface let’s understand what we mean by deformation. There are two broad classes of possible geometric deformations. These go by various names. Uniform deformations (also called affine or linear) includes all modes of deformation in which lines that are parallel prior to the deformation remain parallel after the deformation. There are six types of uniform deformations (Fig. 5).
Figure 5. Alternative modes of uniform shape deformation.
Among these you’ll recognise the deformation modes that are corrected during Procrustes superposition. Nevertheless, the compression/dilation and shear modes can be used, or combined, to describe aspects of genuine shape change. As for the ‘other’ category, it’s usually referred to as a non-uniform deformation in the morphometric literature, but can also be termed a non-affine or non-linear deformation. These are deformations in which lines that are parallel prior to the deformation are not parallel after the deformation. Examples are numerous, but the simplest is the so-called ‘square to kite’ deformation (Fig. 6).
Figure 6. One simple example of a non-uniform shape deformation.
Most ‘real-world’ geometric deformations are combinations of uniform and non-uniform deformation modes. This is certainly the case for the deformation implied by the Acaste and Calymene shape coordinate sets shown in figures 3 and 4.
With our goal of representing shape transformation as a thin plate spline in mind, and with an appreciation of the fact that this spline is (likely) going to be composed of both uniform and nonuniform deformation modes, we’re now in a position to begin a (very generalized) discussion of TPS mathematics. Since we’re going to be minimizing the hypothetical bending energy in the specification of our shape transformation surface, we’re going to need to calculate an index of bending energy at each landmark location. This first step toward this is achieved by the following equation.
(19.1)
Here the value r2ij is the square of the distance between landmarks i and j in the set of shape coordinates for the reference configuration and ln is the natural logarithm function (base e). This calculation quantifies the relative amount of ‘energy’ required to achieve a bend between all pairs of landmarks. The spacing of landmarks represents an important constraint on the spline because it is more difficult (= requires more energy) to achieve a bend between closely spaced landmarks than between landmarks located at a distance from one another. This distinction will have important implications when we discuss principal warps a bit later in this essay series.
Returning to the problem of determining the TPS surface, the various possible modes of bending across a set of landmarks are specified by a partitioned matrix (L) whose structure is summarized as follows.
(19.2)
The P matrix partition summarizes the distances between landmarks using the U function (equation 19.1).
(19.3)
Note here the subscript p refers to the total number of landmarks specified in the reference (and target) configuration(s). The diagonal of this matrix is occupied by zeros because the distance between any landmark and itself is zero. The off-diagonal U values are calculated using equation 19.1. Note that the P the matrix is both square and symmetrical about its diagonal.
The Q matrix summarizes the coordinates of the reference landmark configuration.
(19.4)
The Qt matrix is the transpose of the Q matrix. Finally, the 0 matrix is a 3 x 3 matrix of zeros.
(19.5)
Arranging these matrices in the manner indicated by equation 19.2 produces a composite p+3 x p+3 square, symmetrical matrix whose diagonal elements are all zeros. Once this matrix has been assembled its inverse is calculated (L-1). The L-1 matrix represents the energy required to achieve a displacement of the landmarks of the reference configuration (= bending of the implied reference TPS surface) in any combination and by any amount.
The overall TPS surface is calculated from the L-1 matrix by padding the matrix of x and y shape coordinate values for the target configuration (Xt) out with a 3 x 2 matrix of zeros (Xt+) to give it the same number of rows as the L-1 matrix and then multiplying these two matrices together as follows.
(19.6)
This creates the weight matrix (W). We’ll use the W matrix more in our discussion of partial warps and relative warps. But for now all we need to know is that the W matrix can be partitioned into two sections. The first p rows represent the weights assigned to the non-uniform modes of shape variation and the remaining three rows represent weights assigned to the uniform modes of shape variation.
In order to calculate the 2D representation of the 3D TPS surface the following equations are used.
(19.7) 
(19.8)
The values input into these equations (x,y) are the coordinate positions of vertices of grid centred over the target shape. The dimensions of this grid, the number of grid cells used in its construction, the amount by which the grid extends around the periphery of the target shape landmarks all are under the analyst’s control. Once the positions of the TPS vertices have been established, the grid is constructed by drawing straight lines between adjacent grid vertices to create a deformed, rectilinear grid pattern. By convention, the shape coordinate values of the target shape landmarks are usually plotted along with the grid graphic itself as an aid to interpretation of the TPS surface.
The bending energy matrix is the upper-left p x p block of L-1 matrix (= Lp-1). This partition is used to calculate the non-uniform component of the transformation between the reference configuration and target configuration in the Procrustes space. The equation that yields the coefficients of the non-uniform aspect of the TPS surface is simply the p x 2 matrix product of the target object’s shape-displacement values (Xc) and the p x p bending energy matrix (Lp-1).
(19.9)
This manner of calculating the non-uniform aspect of shape transformation was originally published by Bookstein (1989) and has remained stable since then. Unfortunately, the uniform component of the shape transformation has had a more tortured history. It is beyond the scope of an introductory essay such as this to review that history in detail. The method for estimating the uniform component of shape transformation I will present here is taken from Rohlf (1993) because it is simple and has been used most often in TPS calculations to date. However, there 7 are other methods (see Bookstein 1991, Bookstein 1996, Rohlf and Bookstein 2003, Zelditch et al. 2004).
Under Rohlf’s (1993) approach, and drawing on analogy with the calculation of the bending energy matrix, the uniform component of shape transformation is encoded in the upper right p x 3 block of L-1. This block is referred to by the subscript q, again drawing on analogy with the Q matrix in equation 19.4. Thus, the Lq-1 partition of the L-1 matrix is used to calculate the uniform aspect of the TPS as follows.
(19.10)
Now, let’s run the Acaste and Calymene data through these equations and take a look at their TPS surfaces. Figure 7 summarizes the contrast between the Acaste (Fig. 7A) and Calymene (Fig. 7B) specimens in terms of the ten landmarks used to quantify cranidial morphology. Both of the possible TPS representations are shown, Calymene to Acaste (Fig. 7C) and Acaste to Calymene (Fig. 7D). As per morphometric convention, the configuration of the target specimen is plotted along with the spline. Note the geometrically reciprocal character of two splines.
Figure 7. Total TPS spline surfaces for the geometric transformation between Acaste (red) and Calymene (green). (A) Basal (non-deformed) grid for Acaste landmark configuration. (B) Basal (non-deformed) grid for Calymeme landmark configuration. (C) Thin plate spline surface for the Calymeme-Acaste transformation. (D) Thin plate spline surface for the Acaste-Calymeme transformation.
By way of an example, the Acaste-Calymene TPS surface can be further decomposed into uniform and non-uniform deformation modes, as shown in Figure 8.
Figure 8. Deformational modes of the Acaste-Calymeme geometric transformation. (A) Uniform (affine) TPS surface. (B) Non-uniform (non-affine) TPS surface.
Accordingly, the deformation shown in Figure 7D can be described as a combination of a uniform deformation that combines aspects of a clockwise shear, strong anterio-posterior compression, and clockwise rotation (Fig. 8A). To this is added a pronounced non-uniform deformation centred in the glabella involving a strong element of asymmetrical latero-posterior compression with movement of the three anterior landmarks strongly forward relative to the remaining landmarks. This relative movement results in elongation of the glabella and anterior of the cranidium, lateral migration of the eyes, and latero-posterior migration of the intersection between the crandidium’s the posterior lateral projection and the posterior lateral margin of the glabella (Fig, 8B, see MacLeod 2009, Fig. 5 for landmark definitions).
Thin plate splines for the entire trilobite dataset for which these ten landmarks can be located are shown in Figure 9. In these analyses the sample mean shape was used as the reference shape. Also provided is the value of the total bending energy specified by each spline surface. This number is analogous to the total shape variance and can be used to identify the shapes that deviate more (or less) from the reference (= mean) shape than others.
I hope this brief explanation and demonstration of thin plate splines has demystified the topic for you, at least a bit. Thin plate splines are a very attractive way of graphically depicting shape changes and, because of that, they are also very seductive. They should be used more widely than they are, but they need to be used with caution.
Because these splines are depicted as surfaces that encompass the landmarks themselves, the areas between the landmarks, and even areas outside the region covered by the landmark set, there is tendency to make more of the details of the spline’s configuration than is actually warranted. It should be remembered that, except for the areas immediately surrounding the landmark locations, all other aspects of the spline are artificial interpolations. While the geometry of the spine between the landmarks can identify regions of potential interest (see Bookstein’s 2002 method of TPS creases), interpretations involving these inter-landmark regions should be made with caution. Ideally once inter-landmark regions of interest have been identified, landmarks should (if at all possible) be placed at or near their location and the analysis repeated.
Figure 9. Total TPS surfaces for the comparisons between the trilobite sample mean shape (red) and the landmark shape configuration for 18 genera from the trilobite dataset. Numbers in parentheses beside each genus name are the total bending energies associated with that shape configuration relative to the sample mean shape. This value is analogous to the shape variance.
The issue of the robustness of TPS surfaces to changes in the experimental design should also be mentioned. As in any multivariate analysis, any TPS result is only valid for the specimens used to calculate it, the landmark set used to quantify morphological variation, and the reference form used as a basis for the spline. Variation in any one of these parameters will likely result in substantial changes to the spline’s geometry. In other words, the results of a TPS analysis are simply mathematical descriptions of the shape comparison the analyst has chosen to make between two or more forms; nothing more and nothing less. These descriptions are not generalizable in and of themselves. In subsequent columns we’ll learn how to use the TPS approach to make general statements about the character of shape variation in a sample. But even with those techniques the instability built into the TPS method of shape representation needs to be appreciated. In particular, the selection of the reference shape is, in most cases, a critical decision.
Is the thin plate spline the long-sought realization of the Thompsonian transformation grid concept? In some ways it is and in some ways it isn’t. I suspect Thompson himself would have absolutely loved thin plate splines. D’Arcy Thompson was a great believer in the constraints materials and physical processes place on morphological arrangements. The idea that the TPS algorithm involves a metaphorical concept of bending energy which is required to be minimized by the resulting geometry would have spoken to one of his most deeply held beliefs about the organic world. However, no data or morphological patterns have come to light in the 93 years that have elapsed since On Growth and Form’s publication to lend support the idea that evolutionary processes operate in such a way as to minimize physical parameters such as bending energy. To be sure, organic design cannot exceed the performance limits imposed by the materials used to execute the design. This represents an absolute limitation. But evolutionary history abounds with examples of structures that are inefficient from a purely mechanical point of view. The reason for this this is that mechanical design is only one of the parameters evolutionary processes seek to optimize.
On a more mundane, algorithmic level, the TPS approach also exhibits significant differences with the grids drawn by Thompson and his colleagues, most notably in the sense that Thompson’s grids were conceived of as applying to all points on the form, not just those that happen to be located by particular landmarks. While TPS technique can be used to visualize morphological transformations across and entire form, that transformation is controlled entirely by a relatively small number of point locations. Such a severe abstraction of the overall morphological signal stands in contrast to Thompson’s original transformation-grid concept. Rather, that concept was, as we will see later, much closer in spirit and practice to the analysis 10 of a continuous series of point locations that specify the complex, but biologically informationrich geometries of organic forms. Fortunately, the TPS approach can also be applied to these data despite the fact that the application of TPS analysis to these data have been rare so far because alternatives to TPS-based shape modelling have been available.
A final word about computer programmes for implementing a TPS representation of shape difference. The industry standard remains Jim Rohlf’s tpsSplin and tpsRelw packages (http:// life.bio.sunysb.edu/morph/soft-tps.html). Øvind Hammer’s PAST package (http://folk.uio.no/ ohammer/past/) also calculates thin plate splines though his description of the algorithms it employs to do so (Hammer and Harper 2006) appears to differ in many ways from the canonical descriptions provided by Bookstein (1991), Rohlf (1993) and Zelditch et al. (2004). Other packages that can be used to perform TPS analyses include Dave Sheets’ IMP software (http:// www3.canisius.edu/~sheets/morphsoft.html), Paul O’Higgins’ Morphometrika (http:// sites.google.com/site/hymsfme/downloadmorphologica) and Jon Krieger’s Morpho-Tools online morphometrics data analysis tools site (http://www.morpho-tools.net/).
In his original (1917), and in the later, expanded (1942) edition, Thompson’s goal was to “correlate with mathematical statement and physical law certain of the simpler outward phenomena of organic growth and structure or form while all the while regarding the fabric of the organism, ex hypothesi, as a material and mechanical configuration” (Thompson 1917, p. 17). In other words, Thompson sought to “see how, in some cases at least, the forms of living things, and of parts of living things, can be explained by physical considerations and to realize that in general no organic forms exist save such as are in conformity with physical and mathematical laws” (ibid, p. 15). Thompson was, of course, aware of evolutionary theory and agreed that natural selection operated to sort mechanically efficient from inefficient forms in the manner Darwin had suggested. But he bridled at the idea that every aspect of a form is now, and always has been, under direct adaptive scrutiny, preferring to believe that some aspects of form owe their origin to the physical forces with which they must contend.1 Thompson saw organic form as a ‘diagram of forces’ from which inferences can be made regarding the nature of the forces that act upon it now or that have acted upon it in the past. Using this force metaphor, Thompson saw the mathematical comparison of forms as a way of deducing how these fields of forces changed during both ontogenetic and evolutionary history.
Thompson’s proposed method of force-field analysis was to take two simple line drawings of species’ bodies or some corresponding structural element therefrom (e.g., copepod, ungulate cannon bone, leaf). For convenience we’ll label one form as the ‘reference’ and the other as the ‘target’. In order to better visualize the nature of the geometric transformation Thompson superimposed a rectilinear grid on the reference form. He then worked out simple sets of mathematical transformations that would map point location coordinates of the reference into topologically corresponding locations on the target. Applying those same mathematical transformations to the coordinate locations of the grid vertices Thompson obtained a striking image-based summary of the implied geometric transformation (Fig. 1). In this view Thompson perhaps reflects the same level of discomfort with hyperadaptationist arguments criticized more recently by Gould and Lewontin (1979), among others.
It is clear from Thompson’s many statements throughout the last chapter of On Growth and Form that he regarded the mathematical transformation as pertaining to, and being constrained by, all mathematical points comprising the line drawing and that he respected the principle that biological homology pertained to structures, but not necessarily individual point locations on structures. Rather, it was the configuration of the entire ensemble of mathematical points — represented diagrammatically by the superimposed grid — that he looked to in judging whether he had devised biologically reasonable formulae for a particular form transformation. Similarly, it is clear the only purpose served by the mathematical grid was to passively express the overall geometry of the transformation in the manner of a deformed, map-like coordinate system.
As illustrated in Figure 1, Thompson used his approach to provide examples of both linear (uniform) and non-linear (non-uniform) transformation modes. Like the later ‘relative growth’ studies of Otto Snell, Julian Huxley and Georges Teissier (see below), the thing that impressed morphologists about Thompson’s transformation grids was the fact that seeming complex form changes appeared to be able to be described accurately by simple mathematical transformations applied consistently to all point locations over a form. This suggested to many at the time that the underlying principles and/or determinants of morphological change might be simple when expressed in, or studied using, the language of mathematics.
While Thompson’s transformation-grid approach resulted in the creation of compelling diagrams — so much so both his original set of drawings and many subsequent variations of them have been reproduced in countless books on biology and evolution despite the fact that the physicalforce theory these drawings represent is almost never discussed in those same texts — his geometric approach to the analysis of form never caught on during his lifetime. Thompson himself provided some guidance regarding how to operationalize his transformation grids, which he thought of as a visual tool akin to a modern-day spatial morphing algorithm. Those algorithms subdivide an image into a set of points and then smoothly map a subset of these between a reference and target form with their difference displacements informing the displacement of intermediate points via simple linear interpolation. For example, it is this linear interpolation approach to transformation grid analysis that Thompson used to create his morphlike model of the complex geometric transition between Hyracotherium and modern-day Equus (Fig. 2).
Figure 2. Thompson’s transformation grid analysis for the transition from the Eocene Hyracotherium skull form (A) to the modern Equus skull form (B). Note the representation of hypothetical intermediate stages of the transformation via linear interpolation (C). From Thompson 1917.
Although Thompson’s original, interpolation-based approach to the realization of transformation grids was fine if all you wanted to do was map one form into another, it was not well suited to the summarization of geometric information across a larger sample of data. Arguably Huxley (1924) and Teissier (1929, see also Snell 1892 and Huxley 1932) were more successful in developing an analytic approach to the general problem of form variation than was Thompson. But the logistic regression equations used by students of relative growth — or allometry in modern parlance — were not used to create graphic models of form change with anything like the visual impact of Thompson’s grids.
Curiously, this failure to capitalize on the modelling capabilities of regression-based methods when treating morphological data remained in place for a half-century during which time Huxley and Teissier’s bivariate regression-based approach to form analysis was extended to the multivariate case (via PCA, see Jolicouer and Mosimann 1960) and holistic modelling approaches were developed for other aspects of morphological analysis (see Olson and Miller 1958). As we have seen in the last column, all the mathematical machinery for implementing at least some aspects of useful geometric shape modelling was in place by the 1960s. Yet, no new developments in this area took hold until 1980s despite a few attempts to formulate an explicitly Thompsonian modelling approach in the form of morphological trend surfaces (Sneath 1967) and biorthogonal grids (Bookstein 1978).
In retrospect there appear to be two reasons for this. The first was that the largest school of morphometrics (multivariate morphometrics) tended to look to the communities of statisticians and psychometricians for methodological guidance, neither of which were particularly interested in creating morphological models. The second was that, ever since the 1930s, the tradition in bivariate and multivariate morphological studies was to analyze pairs or sets of distances between landmark locations rather than configurations of Cartesian coordinate locations scattered over a sample of forms. Once the power of shape coordinates had been established by Bookstein (1986) and the outlines of shape theory had begun to emerge (see Kendall 1984), the stage was set to renew the search for an analytic method that could combine the intuitive appeal of Thompson’s transformation grids with the equally popular, and far more powerful, tools of multivariate morphometrics. The key insight that allowed this new approach to shape modelling to be realized was specification of a new spatial metaphor for shape similarity.
In traditional multivariate analysis the similarity between two objects is quantified by calculating the distance between them across all variables. This is fine for a quick-and-dirty summary of form differences, but lacks the ability to track varying patterns of size/shape similarity and difference in different regions of the forms. In the days when morphometricians characterized forms using sets of linear distances between landmarks, this distance-based metaphor seemed both natural and practical. After all, distances are simply magnitudes. There is no information about geometry in a list of distance values. Since geometry can’t be reconstructed precisely from a table of distance values there wasn’t any point in worrying about shape models. But with the move to characterizing forms using the coordinate values of the landmarks themselves — and especially the transformation of landmark coordinate values to Procrustes shape coordinate values via standardization for position, scale, and rotation — it became possible to represent the form and shape similarities or differences between any two objects precisely in a manner that retained the fundamental geometry of the landmark configurations.
Figure 3. Stages in the landmark-based comparison between shapes. (A) Cranidial landmark configuration for the trilobite genus Acaste. (B) Cranidial landmark configuration for the trilobite genus Calymene. (C) Procrustes superposition of the Acaste (red) and Calymene (green) shape coordinates. (D) Shape difference vectors between corresponding reference (red) and target (green) shape coordinate locations. See text for discussion.
Imagine two forms defined by corresponding sets of landmark coordinates. As above, we’ll call one the reference form and the other the target (Fig. 3A-B). Procrustes superposition of these landmark sets transforms the forms into shapes and brings corresponding landmarks into positions of maximal correspondence (= minimal sum of squared deviations, Fig. 3C). The differences between the reference and target shapes can be visualized as a set of difference vectors between the reference and target shapes at each coordinate location (Fig. 3D). Now, rather than summing these differences up to produce a statistical estimate of shape-difference as we’d do in traditional multivariate morphometrics, let’s select the reference shape and express the difference between it and the target shape as a set of vectors with the same displacement as in Fig. 3D, but rotated such that the difference vectors are parallel to the z-axis of a three-dimensional coordinate system. In this system the z-axis expresses the shape difference between the reference and target. The resulting figure (Fig. 4A) expresses the shape difference as a set of stalks at each landmark location. Either the reference or the target shape can be used as the basal shape for constructing this sort of shape-difference diagram. Note that the geometry of the landmark configuration is retained in this graphic expression of shape difference and that we can easily identify which regions of the two shapes are more similar to, or more different from, one another.
Figure 4. Shape-difference diagrams between the Acaste (reference) and Calymene (target) cranidial landmark configurations. (A) Rotation of the reference-to-target shape-difference vectors shown in Fig. 3D to the z-axis to form a three-dimensional representation of shape difference. (B) Spline-estimated surface of the shape-displacement data with colors representing regions of differential curvature. (C) Same spline surface as in B, but this time with the surface represented as a deformed rectilinear grid. In all three diagrams the reference configuration was used as the basis for the shape-difference graphic.
The surfaces shown in Figures 4B and 4C are standard parametric cubic splines. These aren’t terribly useful for summarizing shape difference because they specify elastic deformations in which the shape of the underlying mesh, when viewed from the direction of the shape-displacement (z) axis, is held constant. This amounts to an isomorphic projection of the 3D mesh onto the x,y plane. Instead, the surface interpolation method of choice among morphometricians is an advanced type of polyharmonic spline called a thin plate spline (TPS, see Duchon 1976, 1977). The TPS attempts to mimic the behavior of a defect-free, uniform, and infinitely thin metal plate that is bent in the z-direction to conform to the geometry of the shapedisplacement vectors. This metal-plate metaphor is important because, unlike elastic surfaces, metal plates bend in ways that minimize the energy required to achieve the bend in all directions over the entire plate. In applying this physical metaphor to shape analysis, the TPS represents the surface of minimal bending energy implied by the transformation of one shape into another.
Before we get into the equations for calculating a TPS surface let’s understand what we mean by deformation. There are two broad classes of possible geometric deformations. These go by various names. Uniform deformations (also called affine or linear) includes all modes of deformation in which lines that are parallel prior to the deformation remain parallel after the deformation. There are six types of uniform deformations (Fig. 5).
Figure 5. Alternative modes of uniform shape deformation.
Among these you’ll recognise the deformation modes that are corrected during Procrustes superposition. Nevertheless, the compression/dilation and shear modes can be used, or combined, to describe aspects of genuine shape change. As for the ‘other’ category, it’s usually referred to as a non-uniform deformation in the morphometric literature, but can also be termed a non-affine or non-linear deformation. These are deformations in which lines that are parallel prior to the deformation are not parallel after the deformation. Examples are numerous, but the simplest is the so-called ‘square to kite’ deformation (Fig. 6).
Figure 6. One simple example of a non-uniform shape deformation.
Most ‘real-world’ geometric deformations are combinations of uniform and non-uniform deformation modes. This is certainly the case for the deformation implied by the Acaste and Calymene shape coordinate sets shown in figures 3 and 4.
With our goal of representing shape transformation as a thin plate spline in mind, and with an appreciation of the fact that this spline is (likely) going to be composed of both uniform and nonuniform deformation modes, we’re now in a position to begin a (very generalized) discussion of TPS mathematics. Since we’re going to be minimizing the hypothetical bending energy in the specification of our shape transformation surface, we’re going to need to calculate an index of bending energy at each landmark location. This first step toward this is achieved by the following equation.
(19.1)Here the value r2ij is the square of the distance between landmarks i and j in the set of shape coordinates for the reference configuration and ln is the natural logarithm function (base e). This calculation quantifies the relative amount of ‘energy’ required to achieve a bend between all pairs of landmarks. The spacing of landmarks represents an important constraint on the spline because it is more difficult (= requires more energy) to achieve a bend between closely spaced landmarks than between landmarks located at a distance from one another. This distinction will have important implications when we discuss principal warps a bit later in this essay series.
Returning to the problem of determining the TPS surface, the various possible modes of bending across a set of landmarks are specified by a partitioned matrix (L) whose structure is summarized as follows.
(19.2)The P matrix partition summarizes the distances between landmarks using the U function (equation 19.1).
(19.3)Note here the subscript p refers to the total number of landmarks specified in the reference (and target) configuration(s). The diagonal of this matrix is occupied by zeros because the distance between any landmark and itself is zero. The off-diagonal U values are calculated using equation 19.1. Note that the P the matrix is both square and symmetrical about its diagonal.
The Q matrix summarizes the coordinates of the reference landmark configuration.
(19.4)The Qt matrix is the transpose of the Q matrix. Finally, the 0 matrix is a 3 x 3 matrix of zeros.
(19.5)Arranging these matrices in the manner indicated by equation 19.2 produces a composite p+3 x p+3 square, symmetrical matrix whose diagonal elements are all zeros. Once this matrix has been assembled its inverse is calculated (L-1). The L-1 matrix represents the energy required to achieve a displacement of the landmarks of the reference configuration (= bending of the implied reference TPS surface) in any combination and by any amount.
The overall TPS surface is calculated from the L-1 matrix by padding the matrix of x and y shape coordinate values for the target configuration (Xt) out with a 3 x 2 matrix of zeros (Xt+) to give it the same number of rows as the L-1 matrix and then multiplying these two matrices together as follows.
(19.6)This creates the weight matrix (W). We’ll use the W matrix more in our discussion of partial warps and relative warps. But for now all we need to know is that the W matrix can be partitioned into two sections. The first p rows represent the weights assigned to the non-uniform modes of shape variation and the remaining three rows represent weights assigned to the uniform modes of shape variation.
In order to calculate the 2D representation of the 3D TPS surface the following equations are used.
(19.7) (19.8)The values input into these equations (x,y) are the coordinate positions of vertices of grid centred over the target shape. The dimensions of this grid, the number of grid cells used in its construction, the amount by which the grid extends around the periphery of the target shape landmarks all are under the analyst’s control. Once the positions of the TPS vertices have been established, the grid is constructed by drawing straight lines between adjacent grid vertices to create a deformed, rectilinear grid pattern. By convention, the shape coordinate values of the target shape landmarks are usually plotted along with the grid graphic itself as an aid to interpretation of the TPS surface.
The bending energy matrix is the upper-left p x p block of L-1 matrix (= Lp-1). This partition is used to calculate the non-uniform component of the transformation between the reference configuration and target configuration in the Procrustes space. The equation that yields the coefficients of the non-uniform aspect of the TPS surface is simply the p x 2 matrix product of the target object’s shape-displacement values (Xc) and the p x p bending energy matrix (Lp-1).
(19.9)This manner of calculating the non-uniform aspect of shape transformation was originally published by Bookstein (1989) and has remained stable since then. Unfortunately, the uniform component of the shape transformation has had a more tortured history. It is beyond the scope of an introductory essay such as this to review that history in detail. The method for estimating the uniform component of shape transformation I will present here is taken from Rohlf (1993) because it is simple and has been used most often in TPS calculations to date. However, there 7 are other methods (see Bookstein 1991, Bookstein 1996, Rohlf and Bookstein 2003, Zelditch et al. 2004).
Under Rohlf’s (1993) approach, and drawing on analogy with the calculation of the bending energy matrix, the uniform component of shape transformation is encoded in the upper right p x 3 block of L-1. This block is referred to by the subscript q, again drawing on analogy with the Q matrix in equation 19.4. Thus, the Lq-1 partition of the L-1 matrix is used to calculate the uniform aspect of the TPS as follows.
(19.10)Now, let’s run the Acaste and Calymene data through these equations and take a look at their TPS surfaces. Figure 7 summarizes the contrast between the Acaste (Fig. 7A) and Calymene (Fig. 7B) specimens in terms of the ten landmarks used to quantify cranidial morphology. Both of the possible TPS representations are shown, Calymene to Acaste (Fig. 7C) and Acaste to Calymene (Fig. 7D). As per morphometric convention, the configuration of the target specimen is plotted along with the spline. Note the geometrically reciprocal character of two splines.
Figure 7. Total TPS spline surfaces for the geometric transformation between Acaste (red) and Calymene (green). (A) Basal (non-deformed) grid for Acaste landmark configuration. (B) Basal (non-deformed) grid for Calymeme landmark configuration. (C) Thin plate spline surface for the Calymeme-Acaste transformation. (D) Thin plate spline surface for the Acaste-Calymeme transformation.
By way of an example, the Acaste-Calymene TPS surface can be further decomposed into uniform and non-uniform deformation modes, as shown in Figure 8.
Figure 8. Deformational modes of the Acaste-Calymeme geometric transformation. (A) Uniform (affine) TPS surface. (B) Non-uniform (non-affine) TPS surface.
Accordingly, the deformation shown in Figure 7D can be described as a combination of a uniform deformation that combines aspects of a clockwise shear, strong anterio-posterior compression, and clockwise rotation (Fig. 8A). To this is added a pronounced non-uniform deformation centred in the glabella involving a strong element of asymmetrical latero-posterior compression with movement of the three anterior landmarks strongly forward relative to the remaining landmarks. This relative movement results in elongation of the glabella and anterior of the cranidium, lateral migration of the eyes, and latero-posterior migration of the intersection between the crandidium’s the posterior lateral projection and the posterior lateral margin of the glabella (Fig, 8B, see MacLeod 2009, Fig. 5 for landmark definitions).
Thin plate splines for the entire trilobite dataset for which these ten landmarks can be located are shown in Figure 9. In these analyses the sample mean shape was used as the reference shape. Also provided is the value of the total bending energy specified by each spline surface. This number is analogous to the total shape variance and can be used to identify the shapes that deviate more (or less) from the reference (= mean) shape than others.
I hope this brief explanation and demonstration of thin plate splines has demystified the topic for you, at least a bit. Thin plate splines are a very attractive way of graphically depicting shape changes and, because of that, they are also very seductive. They should be used more widely than they are, but they need to be used with caution.
Because these splines are depicted as surfaces that encompass the landmarks themselves, the areas between the landmarks, and even areas outside the region covered by the landmark set, there is tendency to make more of the details of the spline’s configuration than is actually warranted. It should be remembered that, except for the areas immediately surrounding the landmark locations, all other aspects of the spline are artificial interpolations. While the geometry of the spine between the landmarks can identify regions of potential interest (see Bookstein’s 2002 method of TPS creases), interpretations involving these inter-landmark regions should be made with caution. Ideally once inter-landmark regions of interest have been identified, landmarks should (if at all possible) be placed at or near their location and the analysis repeated.
Figure 9. Total TPS surfaces for the comparisons between the trilobite sample mean shape (red) and the landmark shape configuration for 18 genera from the trilobite dataset. Numbers in parentheses beside each genus name are the total bending energies associated with that shape configuration relative to the sample mean shape. This value is analogous to the shape variance.
The issue of the robustness of TPS surfaces to changes in the experimental design should also be mentioned. As in any multivariate analysis, any TPS result is only valid for the specimens used to calculate it, the landmark set used to quantify morphological variation, and the reference form used as a basis for the spline. Variation in any one of these parameters will likely result in substantial changes to the spline’s geometry. In other words, the results of a TPS analysis are simply mathematical descriptions of the shape comparison the analyst has chosen to make between two or more forms; nothing more and nothing less. These descriptions are not generalizable in and of themselves. In subsequent columns we’ll learn how to use the TPS approach to make general statements about the character of shape variation in a sample. But even with those techniques the instability built into the TPS method of shape representation needs to be appreciated. In particular, the selection of the reference shape is, in most cases, a critical decision.
Is the thin plate spline the long-sought realization of the Thompsonian transformation grid concept? In some ways it is and in some ways it isn’t. I suspect Thompson himself would have absolutely loved thin plate splines. D’Arcy Thompson was a great believer in the constraints materials and physical processes place on morphological arrangements. The idea that the TPS algorithm involves a metaphorical concept of bending energy which is required to be minimized by the resulting geometry would have spoken to one of his most deeply held beliefs about the organic world. However, no data or morphological patterns have come to light in the 93 years that have elapsed since On Growth and Form’s publication to lend support the idea that evolutionary processes operate in such a way as to minimize physical parameters such as bending energy. To be sure, organic design cannot exceed the performance limits imposed by the materials used to execute the design. This represents an absolute limitation. But evolutionary history abounds with examples of structures that are inefficient from a purely mechanical point of view. The reason for this this is that mechanical design is only one of the parameters evolutionary processes seek to optimize.
On a more mundane, algorithmic level, the TPS approach also exhibits significant differences with the grids drawn by Thompson and his colleagues, most notably in the sense that Thompson’s grids were conceived of as applying to all points on the form, not just those that happen to be located by particular landmarks. While TPS technique can be used to visualize morphological transformations across and entire form, that transformation is controlled entirely by a relatively small number of point locations. Such a severe abstraction of the overall morphological signal stands in contrast to Thompson’s original transformation-grid concept. Rather, that concept was, as we will see later, much closer in spirit and practice to the analysis 10 of a continuous series of point locations that specify the complex, but biologically informationrich geometries of organic forms. Fortunately, the TPS approach can also be applied to these data despite the fact that the application of TPS analysis to these data have been rare so far because alternatives to TPS-based shape modelling have been available.
A final word about computer programmes for implementing a TPS representation of shape difference. The industry standard remains Jim Rohlf’s tpsSplin and tpsRelw packages (http:// life.bio.sunysb.edu/morph/soft-tps.html). Øvind Hammer’s PAST package (http://folk.uio.no/ ohammer/past/) also calculates thin plate splines though his description of the algorithms it employs to do so (Hammer and Harper 2006) appears to differ in many ways from the canonical descriptions provided by Bookstein (1991), Rohlf (1993) and Zelditch et al. (2004). Other packages that can be used to perform TPS analyses include Dave Sheets’ IMP software (http:// www3.canisius.edu/~sheets/morphsoft.html), Paul O’Higgins’ Morphometrika (http:// sites.google.com/site/hymsfme/downloadmorphologica) and Jon Krieger’s Morpho-Tools online morphometrics data analysis tools site (http://www.morpho-tools.net/).
REFERENCES
BOOKSTEIN, F. L. 1978. The measurement of biological shape and shape change. Berlin, Springer, 191 pp.
BOOKSTEIN, F. L. 1991. Morphometric tools for landmark data: geometry and biology. Cambridge, Cambridge University Press, 435 pp.
BOOKSTEIN, F. L. 1986. Size and shape spaces for landmark data in two dimensions. Statistical Science, 1, 181–242.
BOOKSTEIN, F., L. 1996. Standard formula for the uniform shape component in landmark data. In L. F. Marcus, M. Corti, A. Loy, G. J. P. Naylor, and D. E. Slice, eds. Advances in Morphometrics. Proceedings of the 1993 NATO Advanced Studies Institute on Morphometrics in Il Ciocco, Italy. New York, Plenum Public Corp., 153–168 pp.
BOOKSTEIN, F. L. 2002. Creases as morphometric characters. In N. MacLeod and P. L. Forey, eds. Morphology, shape and phylogeny. London, Taylor & Francis, 139–174 pp.
DUCHON, J. 1976. Interpolation des fonctions de deux variables suivant le principe de la flexion des plaques minces. R.A.I.R.O. Analyze Numérique, 10(12), 5–12.
DUCHON, J. 1977. Splines minimizing rotation-invariant semi-norms in Sobolev spaces. In W. SCHEMPP and K. ZELLER, eds. Constructive Theory of Functions of Several Variables. Springer, Springer Berlin / Heidelberg. 85–100 pp.
GOULD, S. J. and R. C. LEWONTIN. 1979. The spandrels of San Marco and the Panglossian paradigm: a critique of the adaptationist programme. Proceedings of the Royal Society of London, Series B, 205, 581–598.
HAMMER, Ø. and D. HARPER. 2006. Paleontological data analysis. Oxford, Blackwell Publishing, 351 pp.
HUXLEY, J. S. 1924. The variation in the width of the abdomen of immature fiddler crabs considered in relation to its relative growth-rate. American Naturalist, 58, 468–475.
HUXLEY, J. S. 1932. Problems of relative growth. Methuan & Co., London. 302 pp.
JOLICOUER, P. and J. E. MOSIMANN. 1960. size and shape variation in the painted turtle, a principle component analysis. Growth, 24, 339–354.
KENDALL, D. G. 1984. Shape manifolds, procrustean metrics and complex projective spaces. Bulletin of the London Mathematical Society, 16, 81–121. 11
MacLEOD, N., 2009. Who is Procrustes and what has he done with my data? Palaeontological Association Newsletter, 70, 21–36.
OLSON, E. and R. MILLER. 1958. Morphological integration. Chicago, University of Chicago Press, 317 pp.
ROHLF, F. J. 1993. Relative warp analysis and an example of its application to mosquito wings. In L. F. MARCUS, E. BELLO, and A. GARCÍA-VALDECASAS, eds. Contributions to Morphometrics. Museo Nacional de Ciencias Naturales 8, Madrid. 131–160 pp.
ROHLF, F. J. and F. L. BOOKSTEIN, L. 2003. Computing the uniform component of shape variation. Systematic Biology, 53, 66–69.
SNEATH, P. H. A. 1967. Trend surface analysis of transformation grids. Journal of Zoology, 151, 65–122. SNELL, O. 1892. Die Abhängigkeit des Hirngewichts von dem Körpergewicht und den geistigen Fähigkeiten. Archives of Psychiatry 23, 436–446.
TEISSER, G. 1929. La croissane embryonnaire de Chrysaora hysocella (L.). Archives de Zoologie Expertimentale et Generale, 69, 137–178.
THOMPSON, D. W. 1917. On growth and form. Cambridge University Press, Cambridge. 793 pp. THOMPSON, D. W. 1942. On growth and form. Cambridge University Press, Cambridge. 1116 pp.
ZELDITCH, M. L., D. L. SWIDERSKI, H. D. SHEETS, and W. L. FINK. 2004. Geometric morphometrics for biologists: a primer. Amsterdam, Elsevier/Academic Press, 443 pp.
Endnotes
1 In this view Thompson perhaps reflects the same level of discomfort with hyperadaptationist arguments criticized more recently by Gould and Lewontin (1979), among others.