## Abstract

To characterize the diversity of planar shapes in such instances as insect wings and plant leaves, we present a method for the generation of a smooth morphometric mapping between two planar domains which matches a number of homologous points. Our approach tries to balance the competing requirements of a descriptive theory which may not reflect mechanism and a multi-parameter predictive theory that may not be well constrained by experimental data. Specifically, we focus on aspects of shape as characterized by local rotation and shear, quantified using quasi-conformal maps that are defined precisely in terms of these fields. To make our choice optimal, we impose the condition that the maps vary as slowly as possible across the domain, minimizing their integrated squared-gradient. We implement this algorithm numerically using a variational principle that optimizes the coefficients of the quasi-conformal map between the two regions and show results for the recreation of a sample historical grid deformation mapping of D’Arcy Thompson. We also deploy our method to compare a variety of *Drosophila* wing shapes and show that our approach allows us to recover aspects of phylogeny as marked by morphology.

## 1. Introduction

Relative growth processes in living systems lead to changes in size and shape. Morphometry, the measurement of shape, is often the first descriptive step towards explaining the development and evolution of form in animate and inanimate patterns. Morphometrics, the quantitative study of shape differences between homologous objects, analyses this in terms of a mapping that quantifies the relative positioning of specific landmarks, point-like features in the shapes or on the boundaries of these structures. The identification of such a mapping between two homologous shapes has its roots in the seminal work of D’Arcy Thompson [1], who characterized a ‘theory of transformations’ for a class of planar projections of biological shapes in terms of a grid drawn over a particular shape and then deformed to approximately fit a second, as seen in the example shown in figure 1.

Thompson's ambitions for his grid method were twofold. Firstly, he aimed to uncover the physical basis behind the deformation. In his words [1], once the transformation is known, ‘it will be a comparatively easy task … to postulate the direction and magnitude of the force capable of effecting the required transformation’. The second aim was to ensure that the grid deformation chosen was as elementary as possible while still matching homologous landmarks in the grid coordinates. Thompson further said that whatever process caused the deformation, ‘it is essential that our structure vary in its entirety, or at least that “independent variants” should be relatively few’, citing Occam's razor. Although Thompson provided no quantitative way of calculating the transformations [2], early qualitative approaches included studies of amphibian larval development and ceratopsian dinosaur skulls [3]–[5]. Additionally, the images of deformed grids are useful in pointing out where two homologous structures most differ, using ink-dot experiments on growing leaves [6] or in simulations of neurulation [7] provided. Stimulated by D’Arcy Thompson, Medawar [8] proposed that the creation of deformation grids could be standardized by assuming that the *x*- and *y*-coordinates were transformed according to allometric scaling laws. However, this puts these axes in a privileged position; certainly for an arbitrary shape comparison it is impractical to assume that there will be one set of axes that remains orthogonal under the mapping between them. Sneath [9] also developed a method to calculate a mapping based on the Cartesian grid, but in this case the coordinates of the destination grid were polynomial functions of the source coordinates. The coefficients of the polynomials were calculated to match landmarks under the mapping, and correspond to Thompson's ‘independent variants’. This appears to be the first method to compute mappings, but in essence it is a two-dimensional analogue to a polynomial curve-fitting exercise with coefficients that are hard to relate to local physical or biological quantities. One simple way of constructing a mapping based on landmark data is to connect the points to form a triangular tessellation, and to assume that the mapping is affine in each triangle: this is the finite-element scaling method [10]. A drawback of this mapping is that it is not unique, since there is no canonical way to select the triangular tessellation.

Bookstein's first studies [11] set out to calculate deformation mappings that ensured the proper correspondence between the boundaries of the two structures. His initial calculations approximated the shapes by simple polygons and calculated the correspondence using harmonic maps, ignoring internal constraints. In a later study [12], he chose a mapping that minimized the ‘roughness’ of the deformation subject to internal point constraints (equivalent to demanding that the components of the deformation satisfy biharmonic equations with a point force at each internal landmark). In both cases, the boundary correspondence was prescribed explicitly, by interpolating landmark points on the boundary. Eventually, the idea of creating mappings that were restricted to the interior of the shapes was abandoned in favour of the thin plate spline method [13], eliminating the distinction between boundary and interior landmarks, so that the source and destination objects are defined solely by the positions of these landmark points. This is analogous to the problem of finding the shape of a thin elastic plate, infinite in extent, constrained to prescribed elevations at certain landmark points in the plane, and determined by the solution of the biharmonic equation with a point force at each landmark point. To calculate the thin plate spline mapping between the source and destination shapes, this equation is solved twice: for both the *x*- and *y*-components of the deformation. Where landmark data on the boundary is sparse, one may discretize it and use semilandmarks, whose variation is constrained to be tangential to the boundary. This method has become ubiquitous in the field of morphometry, especially because its linearity has led to a fruitful synthesis with widely used statistical methods, such as kriging [14, 15].

A complementary approach borrows ideas from image registration in medical applications such as MRI and CATscans which involves constructing the deformation that would bring one image in registration with another—precisely the original aim of Thompson. Such methods have led to the field of pattern theory [16, 17] and its relatives such as computational anatomy [18] and have become essential tools for analysing images of highly complex organs, primarily motivated by brain imaging and the calculation of variations in brain shape and the physical characterizations of certain brain pathologies.

These mappings inspired by Thompson's original grid deformation have usually concentrated on only one of Thompson's aims, namely the parsimony of the mapping, or at least that it should be as simple as possible. The second aim, that the mapping should lead to insights into the physical origin of the deformation, has been largely neglected, partly because of the assumption that shape change is directly linked solely to genetic variation. While this is certainly true at one level [19], it has also become increasingly clear that the genetic influence on morphogenesis is mediated by a combination of both developmental and physical processes [20–23]. Indeed, rather generally, there is increasing evidence in a variety of systems that during development, genes influence the distribution of biochemical gradients, or morphogenetic fields [24], which could in turn change the rate of relative growth and rheology of a developing tissue. A growth field in linear elastic or viscoelastic materials, for instance, appears as a distributed body force in the equations of equilibrium, providing a possible interpretation for Thompson's use of the word ‘force’ in describing the physical basis of shape change, and this has become an area of increasing study recently [25, 26].

Thus, an ideal for morphometry is to provide a framework within which descriptive mappings that characterize shape may be married to predictive models of the processes that underlie shapes.Q1 Of the methods described earlier, the closest to this philosophy are the physics-based image registration methods [27, 28], containing distributed force fields as hypothesized by Thompson. For instance, in the method of Bajcsy *et al.* [27], the source image is treated as an elastic body which is deformed by means of a distributed body force, chosen so that the result closely matches the destination image. Alternatively, Christensen *et al.* [28] treat the image being deformed as a viscous fluid. Further examples are provided by Holden [29]. Modern-day implementations of such physics-based registration often take the form of Lagrangian-based deformation calculations [18], where the physics defines an energy satisfied by the (time-dependent) deformation. The deformation is found by solving the partial differential equation (PDE) arising from the minimization of this energy (subject to landmark preservation). However, the mappings here are over the whole image in question rather than restricted to the deforming object itself, so the elastic or viscous properties of the deformation cannot be interpreted as belonging to the object. We may then ask if it is possible to balance the two goals set out by D’Arcy Thompson, can we derive a simple mapping that allows us to efficiently compare shapes while also providing a modicum of insight into the underlying physical processes?

Here we present a method to compare two homologous shapes (*Ω*_{0} and *Ω*_{1}); we will limit ourselves to two dimensions in terms of a mapping from *Ω*_{0} to *Ω*_{1}, with landmark points in *Ω*_{0} mapped to the corresponding points in *Ω*_{1} both for mathematical simplicity, and because there are many nice biological planforms where we can still deploy these methods fruitfully. As desiderata, we would like our mapping to allow us to identify not just points in *Ω*_{0} and the corresponding points in *Ω*_{1}, but also interpolate between the landmark points. Furthermore, given an infinitesimally small circle in *Ω*_{0}, its image under the mapping should tell us about the nature of the mapping at that location: in general the image will be an ellipse, with its orientation characterizing the local rotation, its eccentricity measuring the local shear and the change in area relative to the original circle describing the local expansion or contraction of the material at that point. Since the space of dilatations and shears, combined with rotations completely describes the kinematics of growth and deformations, if our mapping reflects these measures, they may be naturally correlated with predictive models of the shape difference between *Ω*_{0} and *Ω*_{1}.

Mathematically, our mapping between the domains *Ω*_{0} and *Ω*_{1} should allow us to fix a number of points along the boundary and in the interior of the domain, as shown in figure 2, and are characterized in terms of the functions (*u*(*x*,*y*),*v*(*x*,*y*))∈*Ω*_{1}, with (*x*,*y*)∈*Ω*_{0}. A predictive theory of morphometrics would require some dynamical laws for these functions, reflecting their material properties, the biochemical gradients that drive growth or shrinkage, coupled with the balance laws for mass and momentum conservation. The genetic or evolutionary influence on the shape changes is manifested through parameters in these equations. It is in the determination of these parameters that Thompson's parsimony requirement should apply (in contrast to the previously cited models, which choose the deformation itself to be as simple as possible). Thus, given a biophysical model, we seek the simplest parameter field that can explain the shape change observed.

It is natural that our mapping between *Ω*_{0} and *Ω*_{1} be a bijection, mapping the boundaries to each other and preserving the positions of homologous landmarks, with some freedom characterized by a dependence on an underlying parameter field, that encapsulates the influence of evolutionary or genetic processes. In a complete description, these parameters and the mapping would follow from genetics, biochemical and biophysical principles. In the more limited setting considered here, our mapping should reflect as much of these desiderata as possible, and our morphospaces should be described by parameter fields that are the simplest ones able to characterize the difference between the shapes.

Quasi-conformal mappings are relatively simple mappings that satisfy the properties above: they are not as rigid as conformal maps or thin plate splines in that they have free parameters, indeed they have an entire function that is unknown, and yet they are relatively simple to compute numerically. The equations governing such mappings incorporate a parameter field which prescribes the orientation and eccentricity of local patches characterized in terms of ellipses, and thus reflect the natural physical processes associated with planar tissue shaping, but are not as complex as a complete description of the biophysics would entail in terms of material rheology and relative growth rates as embedded in a continuum mechanical field theory. Given the freedom associated with the infinite choices of the parameter field that can describe any given shape change, we postulate our choice of one based on a simple criterion: choose one that varies spatially as little as possible, postulated as an optimization principle.

In §2, we set the stage by describing the basic properties of quasi-conformal maps that we require. In §§3 and 4, we describe the construction of optimal quasi-conformal maps, while in §5, we discuss a numerical method for their determination and present a number of examples that illuminate our theory, drawn from D’Arcy Thompson's work as well as recent work on insect-wing shape.

## 2. A primer on quasi-conformal maps

We start with a brief summary of some of the properties of quasi-conformal maps that we need for our exposition; a complete introduction can be found in the monograph of Ahlfors [30]. Quasi-conformal maps were originally developed as a generalization of the well-known conformal maps by relaxing the condition of conformality that does not allow for any shear deformations, and yet still well suited to treatment by complex analysis. For a general mapping *u*=*u*(*x*,*y*), *v*=*v*(*x*,*y*) between two planar regions, we may immediately write *w*=*u*+*iv* and *z*=*x*+*iy*, so that the mapping can be represented by where . For infinitesimally small regions of (*x*,*y*) space, characterized by the differentials (*dx*,*dy*), the mapping becomes
2.1which is an affine map that takes infinitesimally small circles to infinitesimally small ellipses, with subscripts representing partial differentiation, i.e. *a*_{b}=∂*a*/∂*b*. Equivalently, we have , and so if an infinitesimal circle in *Ω*_{0} is represented as *dz*=*dre*^{iθ}, then
2.2where *w*_{z}=|*w*_{z}|*e*^{iα} and . This equation represents an ellipse whose long axis is at an angle (*α*+*β*)/2 to the horizontal, and with a ratio of long axis to short axis given by
2.3We note that *D*_{w}=1 if and only if , which are just the Cauchy–Riemann equations for conformal maps wherein infinitesimally small circles are mapped to infinitesimally small circles. We may write *D*_{w}=(1+|*μ*|)/(1−|*μ*|), where the quantity is often misleadingly known as the complex dilatation or the first complex dilatation (misleading because in common usage, ‘dilatation’ commonly refers to area changes, whereas *μ* only encodes the eccentricity and orientation of image ellipses), although we will persist with this usage. For later use, we also note that the second complex dilatation is another important parameter, given by , with |*μ*|=|*ν*|.

With these properties defined, a mapping is said to be quasi-conformal if *D*_{w} is bounded, and specifically it is called *K*-quasi-conformal if *D*_{w}≤*K*. It follows immediately that a 1-quasi-conformal map is a conformal map. This condition on *D*_{w} naturally induces a condition on |*μ*| (and on |*ν*|): we need |*μ*|=|*ν*|≤*M*<1 for the map to be defined as quasi-conformal, where *M*=(*K*−1)/(*K*+1). Two further properties of quasi-conformal maps that we will find useful are (i) they can be composed to produce further quasi-conformal mappings; thus the composition of a *K*_{1}-quasi-conformal and a *K*_{2}-quasi-conformal mapping is *K*_{1}*K*_{2}-quasi-conformal, and (i) the inverse of a quasi-conformal map is also a quasi-conformal map.

From the definitions of *μ* and *ν*, we obtain the following differential equations for *w*:
2.4known as the Beltrami equation and the conjugate Beltrami equation, respectively, whose solutions (with the restriction on |*μ*| or |*ν*|) are quasi-conformal maps. We note that if *μ*=*ν*=0, *D*_{w}=1 and hence the mappings are conformal—and indeed in this limit both the Beltrami equations reduce to the Cauchy–Riemann equation .

While either of the equations in (2.4) can be solved to give a quasi-conformal map with the desired shear properties, the conjugate Beltrami equation may be derived naturally from a simple variational principle (presented in §4). This makes it amenable to solve (2.4)_{2}, given *ν* and the boundary of the image (∂*Ω*_{1}) using a Galerkin or finite-element method, and thus find a quasi-conformal map from *Ω*_{0} to *Ω*_{1}. However, we must note the map will not be unique, as we can always compose the computed map with a conformal map, which has three degrees of freedom, and the composite map will have the same distribution of *ν*. With this caveat, we may write , and thus the required quasi-conformal map solution as
2.5and
2.6

## 3. Optimal quasi-conformal mappings

For quasi-conformal mappings between two homologous shapes that reflect both the descriptive and the predictive properties desired in a morphometric measure, we note that the local shear properties of the mapping are controlled by the parameter *ν*=*κ*+*iλ*. If the landmark points in *Ω*_{0} are and are to be matched to in *Ω*_{1} (for *i*=1,…,*N*), then mathematically we would like to find a distribution of *κ*(*x*,*y*), *λ*(*x*,*y*) for (*x*,*y*)∈*Ω*_{0} such that
3.1and
3.2for *i*=1,…,*N*.

However, since *κ* and *λ* are continuous functions of (*x*,*y*) and the 2*N* constraints are discrete, there will be an infinite number of possible quasi-conformal maps. To choose an optimal mapping with characteristics that are mathematically regular and also biologically relevant still leaves us with too many choices. For example, one such map is that of an extremal quasi-conformal map that minimizes the maximum value of . In the absence of point constraints such an optimization would lead to a conformal map (*ν*=0). With these constraints however, minimizing the maximum value of a function over a domain is numerically problematic. To see why, we note that the solution method proceeds by making trial changes in the values of *κ* and *λ* at discretized points in the domain. For most of these trial changes, is unchanged and so the optimization scheme is unable to find an appropriate direction in which to search for a solution. Another possibility is to consider minimizing , which in the limit would tend towards . However, numerical experimentation indicates that optimal solutions of such a process leads to mappings with |*ν*|=*const*. almost everywhere, except in the neighbourhood of the point constraints, where |*ν*| has localized peaks near the point constraints, which is unsatisfactory given that the choice of landmark points can be somewhat arbitrary. This leads us to consider minimizing the gradient of , so that an optimal solution will have constant |*ν*| across the domain, and thus a constant local shear. This is indeed what we find; however, the angle of shear would still vary greatly across the domain as it is not penalized. To prevent this, we finally settled on a proposal that the gradients of both *κ* and *λ* should be minimized, giving rise to the following objective function:
3.3This suggests the following mathematical procedure: find *κ* and *λ*, defined on *Ω*_{0}, that minimizes and satisfy the 2*N* point constraints (3.1) and (3.2). This requires us to find the mappings *u*, *v* given the distributions of the functions *κ*, *λ*. Numerically, the solution is found in the following way. Given a guess for *κ* and *λ*, the quasi-conformal map is calculated, giving *u* and *v*, and hence the destination of the landmarks under the mapping. Then the optimization procedure calculates a new guess for *κ* and *λ* based on minimizing subject to the mapped landmark points being coincident with their actual positions in the destination domain.

It is possible to replace the hard constraints of (3.1) and (3.2) by soft constraints of the form added to the objective functional . However, this would involve making a choice of the parameter *χ*, which would probably be a problem-dependent action.

## 4. Variational principle for quasi-conformal mappings

While quasi-conformal maps can be defined as solutions to a Beltrami equation, this does not directly lead to a computational method of solution. One way to arrive at such a method is to develop a variational principle, similar to that for conformal mappings [31] that can give rise to the formulae (2.5) and (2.6) in weak form, which are amenable to numerical solution via, for example, the finite element method. However, extending the conformal mapping methods to quasi-conformal methods is not obvious. For example, while Mastin & Thompson [32] minimized the functional , where *α*, *β*, *γ* are related to *μ*, together with boundary conditions on ∂*Ω*_{0} and showed that this was equivalent to solving the Beltrami equations, this method works when the geometry of the original domain *Ω*_{0} is a rectangle, with boundary conditions that were simple to implement. However, it is not immediately clear how their method would extend to more general domains such as ours.

Here, we follow a method inspired by Garabedian [33], which does not suffer from this defect, and follows from the observation that the inverse of the equation (2.4)_{1} follows from a variational principle. Since the composition of two quasi-conformal maps is quasi-conformal, the inverse to the solution of the Beltrami equation is also a quasi-conformal map, given by the solution to a conjugate Beltrami equation. We start with the conjugate Beltrami equation which can be rewritten—using *w*=*u*+*iv*, *z*=*x*+*iy* and *ν*=*κ*+*iλ* — as a system of PDEs:
4.1where
4.2
4.3
and
4.4Next, we consider the following functional:
4.5Then, it follows that among all homeomorphisms (*u*,*v*) that map *Ω*_{0} to *Ω*_{1}, those that minimize are also the solutions to the conjugate Beltrami system (4.1), as justified in appendix A.

In order to obtain an expression for *u* and *v* from (4.5), we must obtain the Euler–Lagrange equations associated with the extremization of (4.5), or an equivalent weak form that is amenable to the finite-element method. Letting , the first variation of is given by
4.6where , is an admissible virtual mapping, so that
4.7and *u*, *v* are found by setting for all admissible virtual mappings , which are not independent at the boundary, as shown in §5.

Since the solution of our primary problem—the minimization of —is coupled to the solution of (i.e. the extremization of ), one may ask whether we should solve a constrained variational problem by minimizing with a Lagrange multiplier *t*. However, this is not optimal for two reasons: (i) minimizing would require adding *u* and *v* to *κ* and *λ* as minimization variables, increasing the search space and hence the computation time enormously, and (ii) for any finite *t*, any solution of the extermination problem does not satisfy exactly so that the calculated values of *κ* and *λ* may not bear any relation to the properties of the mapping (*u*,*v*). Of course, this limitation would become less true as *t* increased, and as we recover the original formulation. Since the singular limit is not ideal for numerical computations, we use a nested optimization approach wherein for each trial value of *κ* and *λ*, is minimized first, and the results then used to extremize .

## 5. Numerical implementation

In order to solve the optimization problem introduced above, we first introduce a triangulation of the domain *Ω*_{0} with *N*_{I} interior points and *N*_{B} boundary points. For the finite-element implementation of the weak form equation , we approximate the mapping functions *u*, *v* by piecewise affine functions. The discretized mapping functions *u*^{h}, *v*^{h} (dependent on the mesh size *h*) are
5.1and
5.2Here *u*_{i}, *v*_{i} are the values of *u* and *v* at the interior nodes of the triangulation, and *ϕ*_{i}(*x*,*y*) are the piecewise affine basis functions, being affine on each triangle, taking the value 1 at node *i* and vanishing on every other node. Since the shape of the boundary of *Ω*_{1} is known in the parametric form *u*_{B}(*θ*), *v*_{B}(*θ*), the mapping at the boundary node *k* is thus given by a specific value *θ*_{k} of the parametrizing variable *θ*, multiplied by piecewise affine basis functions *ϕ*_{Bk}(*x*,*y*).

The resulting approximations *u*^{h}, *v*^{h} are thus piecewise affine functions over *Ω*_{0}, and so the gradients ∇*u*^{h} and ∇*v*^{h} are discontinuous, taking a constant value on each element in the triangulated domain. We assume that *κ* and *λ* are likewise approximated by piecewise affine functions,
5.3In this discretized format, the admissible trial mappings , are obtained by allowing for variations in the 2*N*_{I}+*N*_{B} coefficients *u*_{i}, *v*_{i} and *θ*_{k}
5.4and
5.5Substituting these for , in (4.7), we obtain the three finite-element equations for *u*_{i}, *v*_{i} and *θ*_{k}
5.6
5.7
and
5.8for each *m*=1,…,*N*_{I} and *n*=1,…,*N*_{B}, which are 2*N*_{I}+*N*_{B} equations to solve for the same number of coefficients. Since *κ* and *λ* are affine over each triangle in the discretization, the coefficients *A*, *B* and *C* will be nonlinear functions in the same domain. However, we can approximate *A*, *B* and *C* by affine functions generated from their values at each triangle vertex.

Finally, we note that is replaced by its discretized counterpart 5.9

Thus, the optimization process reduces to minimizing by varying *κ*_{j}, *λ*_{j}, *u*_{i}, *v*_{i} and *θ*_{k}, with constraints , , and
5.10and
5.11

We note that the equations (5.6) and (5.7) given by both and are linear in *u*_{i} and *v*_{i}. Therefore, given a set of values for *κ*_{j}, *λ*_{j} and *θ*_{k}, we can calculate the internal values of *u*_{i} and *v*_{i} by inverting the linear system (5.6) and (5.7) thus allowing us to eliminate *u*_{i} and *v*_{i} from the search space for the minimization. It is worth emphasizing that our point matching constraints associated with the landmarks do not distinguish between those defined in the interior and those defined on the boundary, unlike previous approaches. In particular, of the *N* point constraints embodied in (5.10) and (5.11), there will be *N*_{BF} defined on the boundary that force the value of *θ* to be fixed at those points, and *N*_{IF} defined in the interior. Furthermore, if *θ* is fixed, then the equation (5.8) given by does not hold there since this condition is predicated on *θ* being allowed to vary. Thus we consolidate these conditions: define if *θ*_{n} is not fixed, and if *θ*_{n} is fixed. This condition of course depends on the fixed boundary conditions being nodes of the triangulation, but is not difficult to implement in the mesh generation routine.

Thus the optimization process reduces to the following:
5.125.13where *u*^{h}, *v*^{h} are defined in (5.1) and (5.2), with *u*_{i}, *v*_{i} found by solving the linear system .

The triangulation was generated by DistMesh [34], a mesh-generation code (implemented in Matlab), while the optimization algorithm was programmed through an interface to SNOPT, a general-purpose nonlinear constrained optimization program [35].

The optimization algorithm requires a parametrization *u*_{B}(*θ*), *v*_{B}(*θ*) of the boundary curve. The simplest way to encode this is to identify each point on the curve with the angle made between the *x*-axis and the line from the point to a fixed reference point, as seen in figure 2 (inset). This approach only works for star-shaped domains, but this does not appear to be a significant restriction. For many simple shapes, the functions *u*_{B}(*θ*), *v*_{B}(*θ*) are known analytically. For more general shapes, such as the biological shapes analysed later, we can fit *u*_{B}(*θ*) and *v*_{B}(*θ*) to the boundary curve, assuming that the functions are represented as truncated Fourier series (thus preserving periodicity in *θ*). If the boundaries have well-defined corners, the functions are defined piecewise with different Fourier series in each section. In all cases, the endpoints of the curves are fitted exactly. The details of the numerical implementation are given in appendix B.

## 6. Results and applications

To illustrate the results that the method produces, we will discuss three examples: (i) a test problem to calibrate the method, (ii) a recreation of the Thompsonian mapping displayed in figure 1, and (iii) a comparison of the morphology of wings from a number of *Drosophila* species.

The test case corresponds to mapping a circle with unit radius to an ellipse with semi-major axes in a ratio 1 : 2, generated by the affine mapping *u*=2*x*, *v*=*y*, which has and hence . To ensure that our method returns an approximation to this particular choice of coordinates in the mapping, we need to prescribe some fixed points, and we pick four boundary points and two interior points, as indicated in figure 3. For the domain *Ω*_{0}, we use the triangulation in figure 3*a*, so that the deformation calculated by the mapping should yield the mesh in *Ω*_{1}. The triangles have been shaded (both in *Ω*_{0} and *Ω*_{1}) according to their average value of |*ν*|, which is close to for all triangles—as expected.

As our next example, we apply our method to D’Arcy Thompson's transformation theory in the context of two fish species displayed in figure 1. Using the outline of the entire fish is not practical in our case, as it is not a star-shaped domain. Therefore, we omit the fins to obtain a nearly convex shape, which is amenable to analysis. Our boundary fixed points are taken to be the mouth, the anterior ends of the three fins and the two points at which the tail fin meets the body. Interior fixed points are taken to be the eye centre and the upper attachment point of the side fin. The results of the optimization procedure can be seen in figure 4*a*–*b*, and show that the value of |*ν*| over the domain is largely constant, with a slightly greater value near the tail. A different visualization of the mapping is shown in figure 4*c*–*d*, where we have reproduced the original image from Thompson [1] as displayed in figure 1. Just as Thompson superimposed a grid on his image of *Scarus*, we have superimposed a grid of circles, and follow these circles under our numerically calculated optimal quasi-conformal map on the image of *Pomacanthus*. By tracking the positions of the centres of the circles, and distorting them as if they were infinitesimally small, to indicate the shear inherent in the mapping allows us to compare the results of our mapping with Thompson's deformed grid. While we do not expect the grids to match completely, as Thompson was explicitly looking for simple (affine) deformations, such as those giving rise to circular arcs in this case, we do in fact obtain a qualitatively similar mapping, the main difference being that the images of the vertical lines under the mapping have a smaller curvature in *Pomacanthus* than that conjectured by Thompson.

As our final example, we compare wing specimens from different members of the genus *Drosophila*, a favourite of developmental geneticists, and a nice example of a non-trivial flat shape with many variations induced by genetic mutations. We use the database of the wings created by Edwards *et al.* [36] for *Drosophila* species living on the Hawaiian island chain, which have diversified over many millions of years to fill a number of ecological niches. Two example wings (*D. discreta* and *D. primaeva*) are displayed in figure 5*a*–*b*, where we have indicated the landmark points that we preserve under the mapping from one species to another: the intersection of the veins with each other or the wing edge, which are preserved in the *Drosophila* genus. Since the proximal edge of the wing is rather irregular, extending the shape analysis to this region may cause misleading results. However, as our analysis only works for closed shapes, we must provide a curve which will close the boundary and so choose the curve with , where the coefficients are chosen so that the curve passes through all of the leftmost three red points in figure 5*a*–*b*. In figure 5*c*–*f*, we display the results of performing an analysis to find a mapping between the two species *D. discreta* and *D. primaeva*. We observe that |*ν*| appears to be considerably greater on the proximal side of the wing compared with the distal side, and furthermore, there appears to be some contraction in the mapping near the distal edge of the wing.

The change in shape between the two wings reflects the different evolutionary pressures on the two populations that cause differential rates of growth in different parts of the wing. These results shed light on the changes that have occurred in these growth rates as the two species have diverged. We can extend this result to consider mappings between all possible pairs of wings in a given collection. The results for a subset of eight *Drosophila* wings are displayed in figure 6, with four wings from one phylogenetic grouping (*D. montgomeryi*, *D. digressa*, *D. basisetae* and *D. discreta*), and four from a second grouping (*D. cilifera*, *D. truncipenna*, *D. ornata* and *D. adiastola*). Each row represents a source domain (*Ω*_{0}) and each column a destination (*Ω*_{1}). Displaying results in this manner produces a block structure, with 4×4 blocks on the diagonal representing intra-grouping mappings, and off-diagonal blocks representing cross-grouping mappings.

The wing shading in figure 6 according to the value of |*ν*| shows no obvious pattern. However, one observes that the local variation for comparisons within a particular grouping is smaller than for comparisons between the groupings. Thus, we may consider the scaled discrete *H*_{1} norm of the quasi-conformal mapping as a measure of this local variation or shear gradient (where is given in (5.9)). This leads to the heatmap seen in figure 7, where light colors represent low values of , and thus a high match between the shapes. In this context, a match means that the two shapes can be transformed from one to the other by an affine mapping, i.e. a combination of a constant shear, an expansion, a rotation and a translation. While the correspondence is not exact, there is a distinct separation of the two populations according to their phylogenetic grouping (although *D. cilifera* appears to be less similar to the other members of its grouping), with the mappings between wings of the same grouping being closer to affine than cross-grouping mappings. Moreover, figure 6 shows that, in general, the dissimilarity between the two populations is due to the positioning of the two distal internal landmarks.

## 7. Conclusions

Growing laminae such as wings and leaves often lead to fairly reproducible patterns. To characterize the differences between these planar shapes, we have introduced a new method which provides a mapping between two planar regions, while preserving the positions of homologous points and ensuring that the amount of local shear varies as little as possible.

Our method minimizes gradients in shear but does not penalize extensional deformations, thus lying somewhere between a purely descriptive approach and a completely mechanistic predictive one. The former does not respect any mechanistic constraints and is chosen for mathematical convenience. The latter represents the ideal case where some mechanical and chemical laws are placed on changes in the metric of the mapping based on the mechanical, physical or biological properties of the system being analysed. However, for most systems the underlying physical and regulatory causes of the shape changes are not completely known, and so a completely mechanistic approach is currently out of reach. The advantage of our approach is that it retains the notion that the shape comparison should be based on a constrained mapping using a metric constraint (minimization of ) that is a step towards—and yet vastly simpler than—a complete genetic, biochemical and biophysical mechanistic theory. In particular, the quasi-conformal mapping is not as rigid as a conformal mapping, and allows us to vary an underlying parameter field that mimics the influence of intrinsic physical quantities driving the shape change, and could be useful in morphometric studies, where we have a full description of the shape outline, but only sparse data for interior points. This approach should be contrasted with existing mapping techniques such as the thin-plate spline and finite-element scaling method, which select their mapping based on the simplicity of the deformation itself, and are not focused on the processes driving that deformation. Thus, while the simplicity of the existing methods have advantages for statistical analyses, a potential strength of our method is in building a bridge between the evo-devo approaches [19] that link morphometrics to genetics and the biophysical processes which induce shape change.

By applying this method to study two-dimensional biological structures such as insect wings, we can uncover correlations between shapes in terms of properties of the optimal quasi-conformal maps that link them, and thus connect shape-associated phenotypic data with phylogeny and eventually even the genotype. Simultaneously, it allows us to see if any local features of the numerically generated mapping (characterized by the local ellipticity and rotation) correlate with physically observed characteristics, such as cell movement, shape change or tissue stiffness. Through this route one may begin formulating models that capture the physical mechanisms behind the shape differences. As we learn more about mechanisms, our hope is that the statistical and mechanistic approaches will eventually converge, yielding both a descriptive and a prescriptive quantitative view of natural morphospaces.

## Appendix A. Variational principle for quasi-conformal maps

In this section, we provide a justification for the claim that minimizers in the space of homeomorphisms {(*u*,*v*):*Ω*_{0}→*Ω*_{1}} of as given in (4.5) satisfy the conjugate Beltrami system (4.1), following Dierkes *et al.* [37], pp. 263–275 who show that minimizers of the Dirichlet integral in the space of homeomorphisms are conformal mappings, i.e. they satisfy the Cauchy–Riemann system *v*_{y}=*u*_{x}, *v*_{x}=−*u*_{y}.

Setting(*u*_{1},*u*_{2})=(*u*,*v*) and (*x*_{1},*x*_{2})=(*x*,*y*), the inner variation of a functional
A1is given by
A2where *λ*_{α} is an infinitesimal variation made in the independent variable *x*_{α}, with the summation convention implied, and the notation *u*_{κ,λ}=∂*u*_{κ}/∂*x*_{λ}.

For the functional defined in (4.5), we have for the matrix *A*_{αβ} defined by
A3Then the inner variation becomes
A4The minimizing mapping (*u*,*v*) is that for which for all *λ*_{α}. Following the argument in Dierkes *et al.* [37] for conformal maps, this condition is equivalent to
A5Rewriting these, we obtain
A6and
A7Equation (A7) implies that (*u*_{y},*v*_{y})⋅(*Au*_{x}+*Bv*_{x},*Bu*_{x}+*Cu*_{y})=0, or
A8for some . Substituting for *u*_{y} and *v*_{y} in (A6) gives us *χ*^{2}=1, and the root *χ*=+1 is identified by demanding that the Jacobian of the mapping, *u*_{x}*v*_{y}−*u*_{y}*v*_{x}, is positive. We thus obtain the Beltrami system (4.1) as required.

## Appendix B. Optimization problem solution procedure

1.

(a) Obtain parametric representations of the boundary shape for the source (

*x*_{B}(*θ*),*y*_{B}(*θ*)) and destination (*u*_{B}(*θ*),*v*_{B}(*θ*)) shapes, using image analysis techniques.(b) Find the coordinates of the internal landmarks for the source and destination ,

*i*=1,…,*N*_{IF}.(c) Identify the fixed boundary landmarks by finding the source parameter values and the corresponding destination values ,

*k*=1,…,*N*_{BF}.

2. Use (

*x*_{B}(*θ*),*y*_{B}(*θ*)) to find a triangulation of the source domain*Ω*_{0}, with prescribed as nodes of the mesh.3. For each triangle

*a*in the triangulation, and for*k*=1,2,3, calculate the basis functions .4. Use the basis functions and the values of to find the

*N*_{IF}×(*N*_{I}+*N*_{B}) matrix*M*_{ij}such that the destination of under a discretized mapping*u*_{j},*v*_{j}(where*j*are the nodal points) is .5.Initialize the state variables: let

*κ*_{i},*λ*_{i}be zero at each nodal point*i*; let*θ*_{j}be a linear interpolation of the destination parameter values . Ensure*θ*_{1}is a boundary landmark.6.

**Main solution routine.**SNOPT calculates the optimal*y*_{j}(*j*=1,…,*m*) such that*F*_{1}(*y*_{j}) is minimized subject to for*i*=2,…,*n*.(a)

**Subroutine:**calculate*F*_{i}given*y*_{j}.i Let

*κ*_{i}=*y*_{i}and*λ*_{i}=*y*_{NI+NB+i}for*i*=1,…,*N*_{I}+*N*_{B}; let*θ*_{j}=*y*_{2(NI+NB)+j}for*j*=1,…,*N*_{B}.ii Use

*κ*_{i},*λ*_{i}and to calculate .iii Given

*θ*_{j}, find the values of the mapping on the boundary:*u*_{B}(*θ*_{j}) and*v*_{B}(*θ*_{j}).iv Use these with

*κ*_{i},*λ*_{i}, to calculate the matrix and right-hand side of the linear system (5.6)–(5.7).v Invert the system to give the values of

*u*and*v*at the interior nodes. Combine with the boundary values to give*u*_{i},*v*_{i}defined on each node (*i*=1,…,*N*_{I}+*N*_{B}).vi Use

*u*_{i},*v*_{i},*κ*_{i},*λ*_{i}, to calculate for*n*=1,…,*N*_{B}.vii Let

(b) Let

7. Use the output

*κ*_{i},*λ*_{i},*θ*_{k}from the SNOPT routine to calculate the mapping*u*_{i},*v*_{i}by solving the same linear system as in the subroutine.

- Received November 1, 2012.
- Accepted January 24, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.