Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2936781

Formats

Article sections

- Abstract
- Fourier-based Shape Modeling
- Registration
- An Example Application with Damselflies
- Practical Considerations
- LITERATURE CITED

Authors

Related links

Evolution. Author manuscript; available in PMC 2010 September 10.

Published in final edited form as:

Evolution. 2009 April; 63(4): 1003–1016.

Published online 2009 October 17. doi: 10.1111/j.1558-5646.2008.00557.xPMCID: PMC2936781

NIHMSID: NIHMS228585

Associate Editor: G. Hunt

See other articles in PMC that cite the published article.

Quantifying morphological shape is a fundamental issue in evolutionary biology. Recent technological advances (e.g., confocal microscopy, laser scanning, computer tomography) have made the capture of detailed three-dimensional (3D) morphological structure easy and cost-effective. In this article, we develop a 3D analytic framework (SPHARM—spherical harmonics) for modeling the shapes of complex morphological structures from continuous surface maps that can be produced by these technologies. Because the traditional SPHARM methodology has limitations in several of its processing steps, we present new algorithms for two SPHARM processing steps: spherical parameterization and SPHARM registration. These new algorithms allow for the numerical characterization of a much larger class of 3D models. We demonstrate the effectiveness of the method by applying it to modeling the cerci of *Enallagma* damselflies.

The ability to objectively quantify the phenotype of an organism is fundamental to evolutionary biology. Without this ability, we cannot objectively compare phenotypic attributes across taxa or describe how phenotypes change through time. The revolution in quantitative morphometric analyses over the past 20 years has brought tremendous rigor to the quantification and comparison of many types of morphological structures, particularly those with relatively many homologous landmarks (Bookstein 1991; Rohlf and Marcus 1993; Marcus et al. 1996; Small 1996; Dryden and Mardia 1998; Klingenberg et al. 2002; Richtsmeier et al. 2002; Klingenberg and Monteiro 2005). Until recently, a dearth of landmarks on a structure meant that generating high-quality data in three dimensions was difficult, and as a result analytical techniques to exploit such data were unnecessary.

However, recent advances in imaging, microscopy, laser scanning, and computer tomography (CT) have all made the acquisition of rich three-dimensional (3D) morphological data quick and easy from almost any structure, and so analytical techniques to quantify complex 3D shapes are also now sorely needed. A number of approaches have been developed to analyze 3D data, including but not limited to extensions of geometric morphometrics and various spline methods to 3D landmark positions (Bookstein 1991; Rohlf and Marcus 1993; Davis et al. 1997; Dryden and Mardia 1998; Richtsmeier et al. 2002; Zelditch et al. 2004; Gay-Bellile et al. 2006; Liao et al. 2006), eigenshapes of continuous curves (MacLeod 1999), nonlinear registration and surface-to-surface distances (Kristensen et al. 2006; Nieman et al. 2006; Ólafsdóttir et al. 2007), and fractal analyses (Scott et al. 2005) and geographic information systems (GIS) analyses (Plyusnin et al. 2008) of digital elevation surface maps. Each of these techniques has strengths and weaknesses (Richtsmeier et al. 2002; Polly 2008). Techniques such as geometric morphometrics form a rich mathematical description of shape change, but their application requires that a large set of homologous landmarks must be present on every object included in the analysis. Others (e.g., registration and surface-to-surface distance techniques) can analyze collections of objects with few and incomplete sets of landmarks, but are poor descriptors of the underlying objects.

If a morphological structure can be described by a continuous 3D coverage of points (e.g., the voxels corresponding to a structure from a digital image stack produced by imaging or CT scanning), a number of techniques now exist to model the surface from such data, including spherical harmonics (Ballard and Brown 1982; Brechbuhler et al. 1995), hyperquadrics (Hanson 1988), and superquadrics (Terzopoulos and Metaxas 1991). Spherical harmonics, the extension of Fourier techniques to three dimensions, is particularly well suited to modeling shapes from such data, and is being applied to related problems in computer vision (Ballard and Brown 1982; Brechbuhler et al. 1995), computer graphics (Funkhouser et al. 2003; Bulow 2004; Zhou et al. 2004), medical image analysis (Schudy and Ballard 1979; Gerig et al. 2001a,b; Shen et al 2004), and bioinformatics (Ritchie and Kemp 1999). Spherical harmonics were first used as a type of parametric surface representation to construct a functional basis for star-shaped objects (i.e., the object boundary can be defined as a single-valued radius function in a polar coordinate system) (Schudy and Ballard 1979; Ballard and Brown 1982). An extended method, called SPHARM, was proposed to model a much larger class of objects, including those with protrusions, intrusions, and other arbitrarily shaped but simply connected 3D objects (Brechbuhler et al. 1995). SPHARM shares features with the Elliptic Fourier Descriptor (EFD) for describing two-dimensional (2D) contours (Kuhl and Giardina 1982) and can be thought of as a 3D extension of the EFD method.

Spherical harmonics overcomes many of the problems of other techniques: (1) Only a small set of homologous landmarks is required to perform the analyses. These landmarks are used to register objects relative to one another, but not necessarily in the description of the object’s shape. Thus, spherical harmonics is particularly amenable in cases in which geometric morphometrics is not—namely, objects with continuous surfaces but few homologous landmarks. (2) Spherical harmonics produces an orthogonal basis for mathematically representing the 3D shape of a large class of objects. As such, spherical harmonics produces a phenotypic space in which 3D biological structures of heterogeneous shapes can be ordinated, and thus provides a basis for reconstructing the tempo of evolutionary change among those structures.

In this article, we describe the application of the SPHARM algorithmic framework to 3D morphological structures. To build intuition in the approach, we first describe various ways to apply Fourier methods to quantify 2D contours, highlight some of the basic problems with functionally describing shapes in such a framework, and illustrate how these problems can be alleviated. We then describe the SPHARM framework for 3D objects, including algorithmic additions that increase its robustness and applicability to morphological structures. We demonstrate its application using the male mating structures of *Enallagma* damselflies (Odonata: Coenagrionidae). Specifically, we show how SPHARM produces a high-dimensional phenotypic space to describe and quantify the positions of various shapes relative to one another, to perform multivariate descriptive and hypothesis-testing statistical analyses on shape variation, and to estimate ancestral 3D morphologies and rates of evolutionary shape change.

We describe Fourier-based techniques for modeling 2D and 3D shapes. For the purpose of building intuition, we first describe approaches to model closed 2D contours, which closely follows previous descriptions of these methods but with some exceptions. We then show how these techniques extend to modeling 3D surfaces.

A closed 2D contour can be modeled using standard Fourier-based techniques (e.g., the elliptical Fourier analysis ([Rohlf and Archie 1984]). Consider a very simple 2D closed contour—a square (Fig. 1). This contour can be described by a one-parameter function *r*(θ), which maps the distance from a specified origin to each point on the contour as a function of the angle in the polar coordinate system (Fig. 1). The function *r*(θ) can be expressed in terms of a Fourier series

$$r(\mathrm{\theta})=\frac{{a}_{0}}{2}+{\displaystyle \sum _{n=1}^{\infty}({a}_{n}\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}n\mathrm{\theta}+{b}_{n}\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}n\mathrm{\theta})},$$

(1)

where the Fourier coefficients *a _{n}* and

$${a}_{n}=\frac{1}{\mathrm{\pi}{\displaystyle {\int}_{-\mathrm{\pi}}^{\mathrm{\pi}}r(\mathrm{\theta})\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}n\mathrm{\theta}\phantom{\rule{thinmathspace}{0ex}}d\mathrm{\theta}}}$$

(2)

$${b}_{n}=\frac{1}{\mathrm{\pi}{\displaystyle {\int}_{-\mathrm{\pi}}^{\mathrm{\pi}}r(\mathrm{\theta})\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}n\mathrm{\theta}\phantom{\rule{thinmathspace}{0ex}}d\mathrm{\theta}}}.$$

(3)

A square (left panel) can be described by a radius function *r*(θ) (right panel) that maps the distance from a specified origin as a radius sweeps through the square boundary. A few corresponding points are identified on the contour and on the radius **...**

This decomposition provides a generative basis for the function *r*(θ) in terms of sine and cosine functions of varying amplitudes and frequencies. The radial function *r*(θ), and hence the underlying contour, can also be reconstructed to varying degrees of accuracy from this decomposition by using more coefficients in the reconstruction (Fig. 1). With *n* = 0, the radial function is modeled with a single term (i.e., $r(\mathrm{\theta})=\frac{{a}_{0}}{2}$) and the reconstructed contour is, of course, a circle. As more coefficients are added, the contour begins to take form, and with *n* = 100 (i.e., 2*n* + 1 = 201 terms in equation [1]) the reconstruction is almost perfect (the original radial function was sampled at only 256 points).

This Fourier basis representation is not without its limitations. Shown in Figure 2 are three simple contours that do not lend themselves well to the Fourier decomposition described in equation (1). For each shape, the radius crosses the contour at more than one value for many values of θ, given the positions of the origin. Also for the first and third shapes (i.e., the “C” and “L”), the radius does not cross the contour for some angles. Note also that these problems are not necessarily resolved by simply moving the origin. Locations for the origin can be found for the “T” and “L” contours that do resolve these problems, but no location for the origin exists that alleviates these problems for the “C” shape. Several extended Fourier methods (e.g., Kuhl and Giardina 1982; Foote 1989) have been proposed to solve this problem and the rest of this section extends one of these methods (Kuhl and Giardina 1982) that can naturally be extended to model 3D surfaces.

Sample 2D closed contours that cannot be represented by a radius function based on the given origins. No value of *r*(θ) exists for some θs the “C” shape and “L” shapes, and in all of the contours some θs **...**

The 2D Fourier representation can be expanded to these and more complex shapes by employing two parametric functions *x*(θ) and *y*(θ). Similar to the radial function *r*(θ), these functions parameterize the contour as a function of angle θ. However, in this representation, the angle θ is now arc length along the contour specified relative to an arbitrary origin on the contour. The arc length is normalized so that the total length around the contour is 2π, that is the circumference of a unit radius circle. In fact, we can think of this transformation as mapping the contour onto the unit circle, where such a bijective (i.e., one-to-one) mapping preserves the arc length (Fig. 3). The functions *x*(θ) and *y*(θ) correspond to the *x*- and *y*-Cartesian coordinates of the contour. In this parameterization, the “C” contour can now be mapped by two representative functions of θ (Fig. 3). Each function, *x*(θ) and *y*(θ), can now be expanded in terms of a Fourier Series

$$x(\mathrm{\theta})=\frac{{a}_{0}}{2}+{\displaystyle \sum _{n=1}^{\infty}({a}_{n}\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}n\mathrm{\theta}+{b}_{n}\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}n\mathrm{\theta})}$$

(4)

$$y(\mathrm{\theta})=\frac{{c}_{0}}{2}+{\displaystyle \sum _{n=1}^{\infty}({c}_{n}\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}n\mathrm{\theta}+{d}_{n}\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}n\mathrm{\theta})},$$

(5)

where the Fourier coefficients are computed as described in equations (2) and (3). As before, these Fourier coefficients can be used to represent and reconstruct the underlying contour (Fig. 3).

These methods are easily extended to represent the surfaces of closed 3D objects. Shown in Figure 4 is a 3D surface parameterized by the radial function *r*(θ, ). The variable θ is the polar (colatitudinal) coordinate with [0π), and is the azimuthal (longitudinal) coordinate with [0 2π) (Fig. 4A). The function *r*(θ, ) specifies, as a function of two angles, the distance from a specified origin to each point on the surface. The functional representation of the shape in Figure 4B is shown in Figure 4C. As in the 2D case, the radial function embodies the underlying 3D surface, and can be decomposed in terms of a Fourier series and transform. This radial parameterization, however, suffers from the same limitation as described in the previous section—surfaces with zero and more than one value per radius cannot be defined by *r*(θ, ).

Spherical coordinate system that maps out all points on a surface as a function of two angles is illustrated in (A). A star-shaped 3D closed surface is shown in (B), which can be described by a spherical radius function *r*(θ, ) shown in **...**

An analogous change in representation to that used in 2D can be used to parameterize a 3D surface. In this case, three functions, *x*(θ, ), *y*(θ, ), and *z*(θ, ), are needed. Similar to the arc length parameterization in 2D, the angular parameters θ and ϕ correspond to an equal area on the surface. Specifically, a surface is mapped onto a unit sphere under a bijective mapping that minimizes the area and topology distortion. In principle, this mapping is exactly analogous to mapping a 2D contour to a unit circle, but a number of thorny practical problems prevent this mapping from being completely straightforward. However, once mapped onto the sphere, the Fourier transformation on the sphere is straightforward to compute. In the first step, a mapping from the 3D closed surface to a unit sphere is determined. A satisfactory mapping often requires a minimization over some measure of distortion. Three typical distortions are length, angle and area (Floater and Hormann, 2004). A mapping is isometric (length preserving) if and only if it is conformal (angle preserving) and equiareal (area preserving). Thus, an isometric mapping is an ideal mapping without any length, angle, or area distortion. However, isometric mappings only exist in certain cases: for example, the mapping of a cylinder onto a plane. An equiareal mapping is a reasonable alternative to a conformal mapping because it allows each unit area on the object surface to be treated equally by assigning the same amount of parameter space to it. Because an uncontrolled equiareal mapping often contains excessive and undesired angle or length distortions, an ideal spherical mapping should minimize angle and length distortions.

The CALD spherical parameterization algorithm was developed to perform this mapping (Shen and Makedon 2006). The CALD algorithm starts from an initial parameterization and performs local and global smoothing methods alternately until a solution is approached. For this algorithm, the surface is represented as a triangular mesh consisting of vertices (a dense sampling of points on the surface) and faces (connections between adjacent vertices to form a tessellated mesh of triangles). Mesh smoothing relocates vertices on the sphere to improve the parameter mesh quality without changing its topology. Note that smoothing operates only on the parameter mesh (e.g., the spherical parameter mesh shown in Fig. 5C) and does not affect the object mesh (e.g., the mesh describing the surface of the object in Fig. 5B) and its shape. CALD contains three key components. The initial parameterization step is an extension of Brechbuhler’s (1995) method for triangular meshes. A local smoothing method and a global smoothing method are then applied to iteratively improve the quality of the parameterization. The local smoothing step aims to minimize the area distortion at a local submesh by solving a linear system and also to control its worst length distortion at the same time. The global smoothing step calculates the distribution of the area distortions for all the mesh vertices and tries to equalize them over the whole sphere. The overall CALD algorithm combines the local and global methods together and performs each method alternately until a solution is achieved.

An arbitrarily shaped but simply connected 3D closed surface can be described by three spherical functions *x*(θ, ), *y*(θ, ), and *z*(θ, ) based on an underlying spherical parameterization (i.e., a bijective **...**

This mapping applies to genus-zero surfaces (Weisstein 2008)—surfaces and objects with no "holes" or “handles” (e.g., a sphere is a genus-zero object, whereas a torus and a teacup are both genus-one objects). Shown in Figure 5A is a bowl-shaped surface, in Figure 5B a triangular mesh version of this surface, in Figure 5C its mapping onto a sphere, and the corresponding functions *x*(θ, ), *y*(θ, ), and *z*(θ, ) are given in Figure 5D–F. These three functions can now be expressed in terms of the Fourier spherical harmonic (SPHARM) functions (Brechbuhler et al. 1995)

$${Y}_{l}^{m}(\mathrm{\theta},\mathrm{\varphi})=\sqrt{\frac{2l+1(1-m)!}{4\mathrm{\pi}(l+m)!}}\phantom{\rule{thinmathspace}{0ex}}{p}_{l}^{m}(\text{cos}\phantom{\rule{thinmathspace}{0ex}}\mathrm{\theta}){e}^{\mathit{\text{im}}\mathrm{\varphi}},$$

(6)

where ${P}_{l}^{m}$ (cos θ) are the associated Legendre polynomials defined by the differential equation

$${P}_{l}^{m}(x)=\frac{{(-1)}^{m}}{({2}^{l}l!)}{(1+{x}^{2})}^{\frac{m}{2}}\frac{({d}^{l+m})}{({\mathit{\text{dx}}}^{l+m})}{({x}^{2}-1)}^{l}.$$

(7)

Each function is independently decomposed in terms of the spherical harmonics as

$$x(\mathrm{\theta},\mathrm{\phi})={\displaystyle \sum _{l=0}^{\infty}{\displaystyle \sum _{m=-l}^{l}{c}_{\mathit{\text{lx}}}^{m}{Y}_{l}^{m}(\mathrm{\theta},\mathrm{\phi})}},$$

(8)

$$y(\mathrm{\theta},\mathrm{\phi})={\displaystyle \sum _{l=0}^{\infty}{\displaystyle \sum _{m=-l}^{l}{c}_{\mathit{\text{ly}}}^{m}{Y}_{l}^{m}(\mathrm{\theta},\mathrm{\phi})}},$$

(9)

$$z(\mathrm{\theta},\mathrm{\phi})={\displaystyle \sum _{l=0}^{\infty}{\displaystyle \sum _{m=-l}^{l}{c}_{\mathit{\text{lz}}}^{m}{Y}_{l}^{m}(\mathrm{\theta},\mathrm{\phi})}}.$$

(10)

These terms can be bundled into a single vector-valued function

$$\upsilon (\mathrm{\theta},\mathrm{\phi})={\displaystyle \sum _{l=0}^{\infty}{\displaystyle \sum _{m=-l}^{l}{c}_{l}^{m}{Y}_{l}^{m}(\mathrm{\theta},\mathrm{\phi})}},$$

(11)

where υ(θ, ) = (*x*(θ, ), *y*(θ, ), *z*(θ, ))^{T} and ${c}_{l}^{m}={({c}_{\mathit{\text{lx}}}^{m},{c}_{\mathit{\text{ly}}}^{m},{c}_{\mathit{\text{lz}}}^{m})}^{T}$.

The coefficients ${c}_{l}^{m}$ are determined using standard least-squares estimation, which we describe next using *x*(θ, ) as an example. Our goal is to compute the coefficients ${c}_{\mathit{\text{lx}}}^{m}$ up to a user-specified maximum degree *L*_{max}. Assume that an input spherical function *x*(θ, ) is described by a set of spherical samples (θ_{i}, _{i}) and their function values *x _{i}* =

$$\left(\begin{array}{ccccc}\hfill {y}_{1,1}\hfill & \hfill {y}_{1,2}\hfill & \hfill {y}_{1,3}\hfill & \hfill \dots \hfill & \hfill {y}_{1,k}\hfill \\ \hfill {y}_{2,1}\hfill & \hfill {y}_{2,2}\hfill & \hfill {y}_{2,3}\hfill & \hfill \dots \hfill & \hfill {y}_{2,k}\hfill \\ \hfill \vdots \hfill & \hfill \vdots \hfill & \hfill \vdots \hfill & \hfill \hfill & \hfill \vdots \hfill \\ \hfill {y}_{n,1}\hfill & \hfill {y}_{n,2}\hfill & \hfill {y}_{n,3}\hfill & \hfill \dots \hfill & \hfill {y}_{n,k}\hfill \end{array}\right)\phantom{\rule{thinmathspace}{0ex}}\left(\begin{array}{c}\hfill {a}_{1}\hfill \\ \hfill {a}_{2}\hfill \\ \hfill {a}_{3}\hfill \\ \hfill \vdots \hfill \\ \hfill {a}_{k}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill {x}_{1}\hfill \\ \hfill {x}_{2}\hfill \\ \hfill \vdots \hfill \\ \hfill {x}_{n}\hfill \end{array}\right)\phantom{\rule{thinmathspace}{0ex}},$$

where ${y}_{i,j}={Y}_{l}^{m}({\mathrm{\theta}}_{i},{\mathrm{\phi}}_{i})$, *j* = *l*^{2} + *l* + *m* + 1, and *k* = (*L*_{max} + 1)^{2}. Note that we use an indexing scheme that assigns a unique index *j* to every pair (*l*, *m*). Least square fitting is used to solve the above system for (*a*_{1}, *a*_{2}, …, *a _{k}*)

$$\hat{x}(\mathrm{\theta},\mathrm{\phi})={\displaystyle \sum _{l=0}^{{L}_{\text{max}}}{\displaystyle \sum _{m=-1}^{l}{\hat{c}}_{\mathit{\text{lx}}}^{m}{Y}_{l}^{m}(\mathrm{\theta},\mathrm{\phi})\approx x(\mathrm{\theta},\mathrm{\phi}).}}$$

The more degrees (i.e., larger *L*_{max}) one uses, the more accurate the reconstruction *$\widehat{x}$*(θ, ) is. The same least-squares estimation is applied to *y*(θ, ) and *z*(θ, ), and coefficients ${c}_{\mathit{\text{ly}}}^{m}\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}{c}_{\mathit{\text{lz}}}^{m}$ are determined accordingly. After that, the bundled coefficients ${c}_{l}^{m}$ approximate the full underlying surface and can be used to represent and reconstruct the surface (Fig. 7).

When objects are being compared, their SPHARM models must be placed into a common reference system in order for corresponding coefficients in the two models to be comparable and to allow appropriate morphing between them. The goal of registration is to bring the objects as near as possible into positional and rotational alignment with one another. The registration requires the specification of six or more landmarks on each shape, corresponding to "similar" points that will guide the registration. Let *P* = {*$\stackrel{\u20d7}{p}$*_{1}, *$\stackrel{\u20d7}{p}$*_{2}, …, *$\stackrel{\u20d7}{p}$ _{n}*} and

In the first registration step, the following objective function is minimized to solve for the rotation, *R*, and translation, *T*, that brings the specified landmarks into closest positional agreement

$${{\displaystyle \sum _{i=1}^{n}\u301a\Vert {\overrightarrow{q}}_{i}-R{\overrightarrow{p}}_{i}-T\Vert \u301b}}^{2}.$$

(12)

We employ a quaternion-based algorithm (Besl and McKay 1992) to find the rotation and translation that minimizes the above least-squares error (the quaternion-based approach affords a simpler closed-form solution for the above formulation). This is similar to the Procrustes analyses, given that the landmark set roughly defines the homology between objects. Figure 6A–C shows a sample result of aligning an object to the template, where landmarks are shown as colored dots on surfaces.

Step 1 in SPHARM registration: (A) the template, (B) an object before alignment, (C) the same object aligned to the template in the object space. Step 2 in SPHARM registration: (A) the template, (B) an object aligned to the template in the object space, **...**

The registration quality can be evaluated by the root mean squared distance *RMSD* between the corresponding surface parts. Let *S*_{1} and *S*_{2} be two SPHARM surfaces, where their SPHARM coefficients are formed by ${c}_{1,l}^{m}\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}{c}_{2,l}^{m}$, respectively, for 0 ≤ *l* ≤ *L*_{max} and 1*l* ≤ *m* ≤ *l*. The *RMSD* between *S*_{1} and *S*_{2} is (Gerig et al. 2001)

$$RMSD=\sqrt{1/4\pi \phantom{\rule{thinmathspace}{0ex}}{\displaystyle {\sum}_{l=0}^{{{\displaystyle L}}_{\mathrm{max}}}{\displaystyle {\sum}_{m=-1}^{l}{{\displaystyle \Vert {{\displaystyle c}}_{l,t}^{m}-{{\displaystyle c}}_{2,l}^{m}\Vert}}^{2}.}}}$$

(13)

Although using a set of specified landmarks will align an object to the template in Euclidean space, the *RMSD* between these two parametric surfaces will not be minimized unless the surface homology between these models is optimized. In the second registration step, we establish such an optimized surface homology. The homology between SPHARM models is implied by the underlying parameterization: two points with the same parameter pair (θ, ) on two surfaces are defined to be a corresponding pair. In other words, the underlying parameterization defines the landmark labels if one treats a surface as an infinite number of landmarks. Thus, to create an ideal homology, we can simply rotate the parameter net of an object to best match the template’s parameter net. In Figure 6 (D–F), we superimpose a colored mesh onto a SPHARM reconstruction to show its underlying parameterization. On the mesh, the yellow, red, and blue dots indicate the north pole (0, 0), the south pole (0, π), and the crossing point of the zero meridian and the equator $(0,\frac{\mathrm{\pi}}{2})$, respectively. The template and its parameterization are shown in Figure 6D, and the initial parameterization of the object to be aligned is shown in Figure 6E. Even though the object has been rotated and translated so that its landmark positions best match those of the template, the *RMSD* (=79.0) between these two models is large because their parameterizations are not well aligned. Rotating the parameterization of the object to the position shown in Figure 6F minimizes the *RMSD* (=12.5) between them. The second step of SPHARM registration aims to achieve this minimized *RMSD* by rotating the parameterization of the object surface.

One solution for rotating the parameterization of a SPHARM model is to recalculate the SPHARM coefficients using the rotated parameterization, but this requires solving three linear systems and is time-consuming. To accelerate the process, we use a rotational property in the harmonic theory and rotate SPHARM coefficients without recalculating the SPHARM expansion. Let $\upsilon (\mathrm{\theta},\mathrm{\phi})={\displaystyle {\sum}_{l=0}^{\infty}{\displaystyle {\sum}_{m=-1}^{l}{c}_{l}^{m}{Y}_{l}^{m}(\mathrm{\theta},\mathrm{\phi})}}$ be a SPHARM parametric surface. After rotating the parameter net on the surface in Euler angles (αβγ), the new coefficients ${c}_{l}^{m}$ (αβγ) are (Bruel and Hennocq 1995; Ritchie and Kemp 1999):

$${c}_{l}^{m}(\mathrm{\alpha}\mathrm{\beta}\mathrm{\gamma})={\displaystyle \sum _{n=-1}^{l}{D}_{\mathit{\text{mn}}}^{l}(\mathrm{\alpha}\mathrm{\beta}\mathrm{\gamma}){c}_{l}^{n}},$$

(14)

where

$${D}_{\mathit{\text{mn}}}^{l}(\mathrm{\alpha}\mathrm{\beta}\mathrm{\gamma})={e}^{-i\mathrm{\gamma}n}{d}_{\mathit{\text{mn}}}^{l}(\mathrm{\beta}){e}^{-i\mathrm{\alpha}m}$$

and

$${d}_{\mathit{\text{mn}}}^{l}={\displaystyle \sum _{t=\text{max}({0}_{l}n-m)}^{\text{min}(l+n,l-m)}{(-1)}^{t}\frac{\sqrt{(l+n)!(l-n)!(l+m)!(1-m)!}}{(l+n-t)!(l-m-t)!(t+m-n)!t!}}\phantom{\rule{thinmathspace}{0ex}}\times \phantom{\rule{thinmathspace}{0ex}}{\left(\frac{\text{cos}\phantom{\rule{thinmathspace}{0ex}}\mathrm{\beta}}{2}\right)}^{(2l+n-m-2t)}{\left(\frac{\text{sin}\phantom{\rule{thinmathspace}{0ex}}\mathrm{\beta}}{2}\right)}^{(2l+m-n)}.$$

The goodness of the match is measured by the *RMSD* between two models, which can be calculated directly from SPHARM coefficients according to equation (13).

We employ a sampling-based strategy that fixes one parameterization and rotates the other to optimize the surface correspondence by minimizing the *RMSD* defined in equation (13). The rotation space can be sampled nearly uniformly using icosahedral subdivisions. This assigns rotation angles to β and γ. Let *n* be the number of icosahedral samples. The expansion coefficients are then rotated through (0βγ) and then by *n* equal steps in α using equation (14), evaluating the *RMSD* at each orientation. The result is the best orientation that minimizes the *RMSD*.

We have employed the 3D SPHARM representation described in the previous sections to model the morphology of cerci, the male morphological mating structures used by males to grasp females during courtship and copulation, and the females where these male structures contact the female in the *Enallagma* damselfly clade (McPeek et al. 2008, 2009). *Enallagma* females use the 3D shape of these structures to assess whether males are conspecifics or heterospecifics for mating (Paulson 1974; Robertson and Paterson 1982), and likewise the 3D shapes of these structures are used by humans to identify males to species (Westfall and May 2006). In that study, we modeled the cerci of 41 *Enallagma* species, and used the spherical harmonic representations to reconstruct the evolutionary tempo and mode of shape change in this important morphological structure across the 15 million year history of the clade. Here, we emphasize the application and discuss the issues of applying spherical harmonics to the evolution of 3D morphological structures using *Enallagma* as examples.

The basic data needed to apply these techniques are a dense sampling of points on the surfaces of the structures to be modeled. Many different technologies are now available to obtain such samplings, and a thorough review of them is beyond the scope of this article. Two that we have used are confocal microscopy and CT. We employed micro CT scanning to generate the data for *Enallagma* cerci at a resolution of 2.5 µm using a SkyScan 1172 high-resolution micro-CT scanner (SkyScan, Kontich, Belgium). The output from the CT scanning was a digital image stack of cross-sections of a damselfly male. Amira™ version 4.1.2 software (Mercury Computer Systems Inc., Chelmsford, MA) was used to identify all voxels associated with each cercus (damselflies have two, a left and a right cercus, which are mirror images of one another, to grasp females) to segment it from the rest of the body. A high-resolution triangular mesh (10,002 vertices forming 20,000 triangles) of the cercus surface was then constructed from the object defined by the segmented voxels. The positions of seven landmarks were then manually identified for registration. Because the spherical harmonic coefficients are influenced by both size and shape, we size-standardized all cerci to a common length between two of the landmarks, thereby making shape the primary difference between the cerci of different species.

First, a spherical parameterization must be constructed for each cercus. The triangular mesh and landmarks for a cercus of *E. anna* are shown as the raw mesh and as the mapping to the unit sphere in Figure 7. To illustrate how greater degrees of spherical harmonic representation give more detailed description of the object, we show reconstructions of this *E. anna* cercus using *l* = 1, 8, 16, and 24 in Figure 7.

Registration is the next important step in any analysis in which shapes are to be compared. We recommend that one of the objects be chosen as the registration “template” and all objects be registered to this template before further analyses. For our analyses, we arbitrarily chose *E. anna* as the template. Figure 8 illustrates the registration steps using *E. antennatum* being registered to the template *E. anna*. First, the object is translated and rotated to bring the landmarks into as close a positional correspondence as possible (Fig. 8A–C), and then the parameter alignment is performed (Fig. 8D–F). This registration step completes a generalized Procrustes analysis as well as optimizes the homology among objects. Figure 9 shows the registered model reconstructions for the cerci of several *Enallagma* species.

Step 1 in SPHARM registration of *Enallagma antennatum* to *E. anna* as a template: (A) the *E. anna* template, the *E. antennatum* cercus (B) before alignment, and (C) after aligning to the template landmarks in the object space. Step 2 in SPHARM registration: **...**

Once all the models have been aligned, the SPHARM coefficients are an orthogonal basis for a high-dimensional phenotype space that ordinates each of the objects relative to one another. If the objects have been standardized to a common size before the SPHARM models are constructed, this ordination primarily describes shape differences among the objects. The registration makes the coefficients commensurable, and thus the coefficients are quantitative representations of the shapes of the objects that can be used to calculate (1) the distance between the shapes of two structures, (2) the angles between the shapes of two objects relative to a third reference object, and (3) moments of distributions for the collection of shapes in a dataset (e.g., the mean and variance of shape). Standard multivariate statistical techniques (e.g., Morrison 2004) can then be applied to describe and test hypotheses about the collection of objects in the dataset. For example, principal components analysis can be applied to reduce the dimensionality of data, describe the major axes of shape variation in this high-dimensional space, and visualize the positions relative to one another (see McPeek et al. (2008, 2009) for an application to the *Enallagma* cerci). For example, the first four principal components extracted from the 41 *Enallagma* cerci accounted for 75.6% of the total variation in shape (McPeek et al. 2008). Figure 10 presents visual representations of the shape variation along these first four principal components.

Visualization of the shape variation along the first four principal components (PCs). Each row shows eigenshapes (i.e., reconstructions of eigenmodes) along a principal axis space along a range from −2 SD (−2 × SD) to +2 SD (2 **...**

The SPHARM models also can be morphed from one to another, because the homology between the two has been established by their optimally aligned underlying parameterizations. Given two SPHARM descriptors *A* and *B*, we can calculate a new descriptor *C* = (1 − ω) *A* + *w B* : *C* is simply a weighted average of *A* and *B*. As ω is changed from 0 to 1, *C* will transfer from *A* to *B* accordingly and form a morphing result between *A* and *B*. In Figure 11 we show such a morphing example from *E. anna* to *E. antennatum*. In an evolutionary context, morphing is directly applicable in ancestral reconstruction techniques—ancestral reconstruction is simply morphing according to an evolutionary model of character change. Various techniques are now used widely to reconstruct the ancestral character states of quantitative phenotypes over the phylogeny of a clade (e.g., squared change parsimony, minimum evolution, evolutionary contrasts analyses: see review in Rohlf (2001)). All these techniques estimate the phenotypes associated with nodes in a phylogeny by calculating some weighted average of surrounding nodes in the tree. All of these weighted averages are simply different specifications of the above morphing function. Morphing can potentially be used to predict the shape of a common ancestor for a set of species (Rohlf 2001) and thus to help reconstruct the evolutionary tempo of shape changes (see McPeek et al. 2008, 2009 for the application of evolutionary contrasts analyses and character reconstructions to the *Enallagma* clade).

Three-dimensional morphological shape has been one of the most elusive aspects of an organism’s phenotype to accurately and quantitatively characterize in all its subtlety. Fourier methods are a venerable mathematical toolkit for such problems and as such have been used for some time to characterize shape in 2D (Rohlf and Archie 1984; Ferson et al. 1985; Foote 1989). Fourier representations do, however, have their limitations. A major criticism of Fourier-based methods is that coefficients do not represent homologous, or in fact any specific, feature of a shape (Bookstein et al. 1982; Ehrlich et al. 1983). Any particular localized change in a shape (e.g., a local protrusion or invagination) will cause simultaneous changes in many Fourier coefficients. As a result, Fourier coefficients cannot be used as characters in a phylogenetic analysis (Zelditch et al. 1995). However, the inability to use them as phylogenetic characters does not negate the utility of Fourier coefficients as effective descriptors of shape, as an effective basis for comparing differences in shape among a collection of objects (MacLeod 1999; Polly 2008), or as a basis for reconstructing the evolution of shape change (McPeek et al 2008, 2009).

Polly (2008) has lucidly discussed many of the theoretical and practical problems for understanding the evolution of structures that are characterized by such homology-free representations as Fourier coefficients. As he illustrates, direct linear changes in genetics and the underlying developmental bases for structures will not necessarily cause linear changes between structures in our mathematical representations of those structures (Polly 2008). This is an argument about how one structure is morphed into another. In fact, in our approach, the registration algorithm implicitly establishes homology between surfaces by optimally aligning their underlying parameterizations, which is similar to assigning appropriate landmark labels to surface samples. This is why the coefficients of two shapes are commensurable. However, the spherical harmonic coefficients do not represent isolated features on the structure, and so individual coefficients cannot be associated with specific features of the structure.

Alignment of the parameter meshes does, however, generate correspondence of spatial locations across structures. Reconstructing shapes from the spherical harmonic coefficients is done on some distribution of vertices across the surface of the unit sphere (e.g., Fig. 9). Alignment of the parameter meshes of two objects means in practical terms that each vertex on one reconstructed shape corresponds to a vertex on the other object (e.g., Fig. 8). Thus, one could compare shapes by reconstructing the (*x,y,z*) coordinates of each shape from the aligned parameter meshes on a common distribution of vertices on the unit sphere. Each vertex is then a corresponding position on the surfaces of all objects—although we stress that these are not the same as homologous landmarks in geometric morphometrics, because landmarks defined a priori on the surfaces of objects (e.g., our registration landmarks) will not be identical vertices on all objects, although they will be very close to one another. Because the mapping from the spherical harmonic coefficients to spatial information is linear, these two representations of the shapes (i.e., spherical harmonic coefficients vs. spatial locations of vertices) contain identical information.

To show this, we reconstructed spatial models of the 41 *Enallagma* cerci using an icosahedron subdivided to level 4 on the unit sphere (giving 5120 triangles [Teanby 2006]), and extracted principal components from the resulting coordinates of these reconstructions. The eigensystem of this principal components analysis was identical to the eigensystem extracted from the spherical harmonic coefficients used to make the reconstructions (i.e., corresponding eigenvalues accounted for the same amount of variation, and principal component scores from the two analyses were [within the bounds of rounding error] perfectly correlated). Thus, the spherical harmonic coefficients can be used to construct surfaces, and the distances between the locations of each vertex of reconstructed shapes can be used to identify localized areas that differ to greater or lesser extents. This is comparable to registration and surface-to-surface distances methods (e.g., Kristensen et al. 2006; Nieman et al. 2006; Ólafsdóttir et al. 2007) but is done on the reconstructed objects and not on the original data. Parameterizing objects using larger values of *L*_{max} would reduce the difference between the original data and the reconstructed object, but in practice even moderate values of *L*_{max} (we used *L*_{max} = 18 for this analysis) should provide adequate descriptors of shapes for this purpose. In addition, spherical harmonics provides the added benefit of quantitatively describing the shape of each object. We are continuing to develop visualization tools to aid these analyses, methods that provide more exact correspondence between vertices on different objects, and links between various types of analyses (e.g., combining SPHARM with 3D implementations of thin plate splines).

We have also described a very simple morphing function (i.e., a linear averaging) that is equally applied to all Fourier coefficients (e.g., McPeek et al. 2008, 2009). One can imagine different morphing functions being applied to different coefficients (e.g., altering the high and low frequency coefficients according to different functions). At this point, we do not know what the correct morphing function is. Moreover, the correct function should reflect the underlying genetic and developmental mechanisms that define the structure’s shape (Steppan 2002; Rice 2004; Polly 2008). However, we have such information for only a small handful of simple traits, let alone complex shapes, and yet evolutionary morphing (e.g., evolutionary contrasts analyses and ancestral character reconstruction: Felsenstein 1985; Rohlf 2001) is a powerful tool in comparative biology.

Recent technological advances in imaging and CT have made capturing the detailed 3D structure of morphological features as data quick and inexpensive. The extension of Fourier methods to such 3-D data is relatively straightforward, and the algorithms we describe here alleviate many of the computational problems inherent in this transition. Further work is obviously needed to continue to refine the SPHARM procedures, but these refinements will continue to develop familiar and powerful analytical tools to capture the rich structure inherent in 3D morphologies.

We would like to thank three anonymous reviewers for comments on previous drafts of the article. This work was supported by NIH/NIBIB R03EB008674-01 to L. Shen and by National Science Foundation grant IBN-0516104 to M. A. McPeek and H. Farid. Software coded in Matlab (The Mathworks, Inc., Natick, MA) is available from the authors on request (E-mail: ude.htuomtrad@keepcm.kram). Models of *Enallagma* cerci are available at http://www.enallagma.com/cerci/damselflyMating.html.

- Ballard DH, Brown CM. Computer vision. Englewood Cliffs, NJ: Prentice-Hall; 1982.
- Besl PJ, McKay ND. A method for registration of 3-D shapes. IEEE Trans. Patt. Anal. Mach. Intelligence. 1992;14:239–256.
- Bookstein FL. Morphometric tools for landmark data : geometry and biology. Cambridge, England, New York: Cambridge Univ. Press; 1991.
- Bookstein FL, Strauss RE, Humphries JM, Chernoff B, Elder RL, Smith GR. A comment upon the uses of Fourier methods in systematics. Syst. Zool. 1982;31:85–92.
- Brechbuhler C, Gerig G, Kubler O. Parameterization of closed surfaces for 3D shape description. Comp. Vis. Image Understanding. 1995;61:154–170.
- Bulow T. Spherical diffusion for 3D surface smoothing. IEEE Trans. PAMI. 2004;26:1650–1654. [PubMed]
- Burel G, Hennocq H. Determination of the orientation of 3D objects using spherical harmonics. Graph. Models Image Proc. 1995;57:400–408.
- Davis MW, Flamig DP, Harms SE. A physics-based coordinate transformation for 3-D image matching. IEEE Trans. Med. Imag. 1997;16:317–328. [PubMed]
- Dryden IL, Mardia KV. Statistical shape analysis. New York: John Wiley & Sons; 1998.
- Ehrlich R, Pharr RB, Jr., Healy-Williams N. Comments on the validity of Fourier descriptors in systematics: a reply to Bookstein et al. Syst. Zool. 1983;32:202–206.
- Felsenstein J. Phylogenies and the comparative method. Am. Nat. 1985;126:1–25.
- Ferson S, Rohlf FJ, Koehn RK. Measuring shape variation of two-dimensional outlines. Syst. Zool. 1985;34:59–68.
- Floater MS, Hormann K. Multiresolution in geometric modelling. Springer; 2004. Surface parameterization: a tutorial and survey.
- Foote M. Perimeter-based Fourier analysis: a new morphometric method applied to the
*Trilobite cranidium*. J. Paleontol. 1989;63:880–885. - Funkhouser T, Min P, Kazhdan M, Chen J, Halderman A, Dobkin D, Jacobs D. A search engine for 3D models. ACM Trans. Graphics. 2003;22:83–105.
- Gay-Bellile V, Perriollat M, Bartoli A, Sayd P. Image registration by combining thin-plate splines with a 3D morphable model. IEEE Image Proc. 2006:1069–1072.
- Gerig G, Styner M, Jones D, Weinberger D, Lieberman J. Shape analysis of brain ventricles using SPHARM. IEEE Workshop on Mathematical Methods in Biomedical Image Analysis; 2001a. pp. 171–178.
- Gerig G, Styner M, Shenton ME, Lieberman JA. Shape versus size: improved understanding of the morphology of brain structures. 4th International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI 2001); Ultrecht, the Netherlands: Springer; 2001b. pp. 24–32.
- Hanson AJ. Hyperquadrics: smoothly deformable shapes with convex polyhedral bounds. Comput. Vision Graph. Image Proc. 1988;44:191–210.
- Klingenberg CP, Monteiro LR. Distances and directions in multidimensional shape spaces: implications for morphometric applications. Syst. Biol. 2005;54:678–688. [PubMed]
- Klingenberg CP, Barluenga M, Meyer A. Shape analysis of symmetric structures: quantifying variation among individuals and asymmetry. Evolution. 2002;56:1909–1920. [PubMed]
- Kuhl FP, Giardina CR. Elliptic Fourier features of a closed contour. Comput. Graph. Image Proc. 1982;18:236–258.
- Liao S-H, Tong R-F, Dong J-X. 3D whole tooth model from CT volume using thin-plate splines. IEEE Computer Supported Cooperative Work in Design. 2005:600–604.
- MacLeod N. Generalizing and extending the eigenshape method of shape space visualization and analysis. Paleobiology. 1999;25:107–138.
- Marcus LF, Corti M, Loy A, Naylor GJP, Slice DE, editors. Advances in morphometrics. New York: Plenum; 1996.
- McPeek MA, Shen L, Torrey JZ, Farid H. The tempo and mode of 3-dimensional morphological evolution in male reproductive structures. Am. Nat. 2008;171:E158–E178. [PubMed]
- McPeek MA, Shen L, Farid H. The correlated evolution of 3-dimensional reproductive structures between male and female damselflies. Evolution. 2009 in press. [PubMed]
- Morrison DF. Multivariate statistical methods. New York: McGraw-Hill; 2004.
- Nieman BJ, Flenniken AM, Adamson SL, Henkelman RM, Sled JG. Anatomical phenotyping in the brain and skull of a mutant mouse by magnetic resonance imaging and computed tomography. Physiol. Genom. 2006;24:154–162. [PubMed]
- Ólafsdóttir H, Darvann TA, Hermann NV, Oubel E, Ersbøll BK, Frangi AF, Larsen P, Perlyn CA, Morriss-Kay GM, Kreiborg S. Computational mouse atlases and their application to automatic assessment of craniofacial dysmorphology caused by the Crouzon mutation
*Fgfr2*. J. Anat. 2007;211:37–52. [PubMed]^{C342Y} - Paulson DR. Reproductive isolation in damselflies. Syst. Zool. 1974;23:40–49.
- Plyusnin I, Evans AR, Karme A, Gionis A, Jernvall J. Automated 3D phenotype analysis using data mining. PLoS One. 2008;3:e1742. [PMC free article] [PubMed]
- Polly PD. Developmental dynamics and G-matrices: can morphometric spaces be used to model phenotypic evolution? Evol. Biol. 2008;35:83–96.
- Rice SH. Developmental associations between traits: covariance and beyond. Genetics. 2004;166 [PubMed]
- Richtsmeier JT, DeLeon VB, Lele SR. The promise of geometric morphometrics. Yearbook Phys. Anthropol. 2002;45:63–91. [PubMed]
- Ritchie D, Kemp G. Fast computation, rotation, and comparison of low resolution spherical harmonic molecular surfaces. J. Comp. Chem. 1999;20:383–395.
- Robertson HM, Paterson HEH. Mate recognition and mechanical isolation in Enallagma damselflies (Odonata: Coenagrionidae) Evolution. 1982;36:243–250.
- Rohlf FJ. Comparative methods for the analysis of continuous variables: geometric interpretations. Evolution. 2001;55:2143–2160. [PubMed]
- Rohlf FJ, Archie JW. A comparison of Fourier methods for the description of wing shape in mosquitoes (Diptera: Culicidae) Syst. Zool. 1984;33:302–317.
- Rohlf FJ, Marcus L. A revolution in morphometrics. Trends Ecol. Evol. 1993;8:128–132.
- Schudy R, Ballard D. Towards an anatomical model of heart motion as seen in 4-D cardiac ultrasound data. The 6th Conference on Computer Applications in Radiology and Computer-Aided Analysis of Radiological Images.1979.
- Scott RS, Ungar PS, Bergstrom TS, Brown CA, Grine FE, Teaford MF, Walker A. Dental microwear texture analysis shows within-species diet variability in fossil homonins. Nature. 2005;436:693–695. [PubMed]
- Shen L, Makedon F. Spherical mapping for processing of 3-D closed surfaces. Image Vision Comput. 2006;24:743–761.
- Shen L, Ford J, Makedon F, Saykin A. A Surface-based approach for classification of 3D neuroanatomic structures. Intell. Data Anal. 2004;8:519–542.
- Small CG. The statistical theory of shape. New York: Springer; 1996.
- Steppan SJ, Phillips PC, Houle D. Comparative quantitative genetics: evolution of the G matrix. Trends Ecol. Evol. 2002;17:320–327.
- Teanby NA. An icosahedron-based method for even binning of globally distributed remote sensing data. Comput. Geosci. 2006;32:1442–1450.
- Terzopoulos D, Metaxas D. Dynamic 3D models with local and global deformations: deformable superquadrics. IEEE Trans. PAMI. 1991;13:703–714.
- Weisstein EW. Genus. From MathWorld A Wolfram Web Resource. 2008. http://mathworld.wolfram.com/Genus.html.
- Westfall MJ, May ML. Damselflies of North America. 2nd edn. Gainesville, FL: Scientific Publishers; 2006.
- Zelditch ML, Fink WL, Swiderski DL. Morphometrics, homology, and phylogenetics: quantified characters as synapomorphies. Syst. Biol. 1995:179–189.
- Zhou K, Bao H, Shi J. 3D surface filtering using spherical harmonics. CAD. 2004;36:363–375.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |