(PDF) Evolutionary morphing

Evolutionary Morphing David F. Wiley∗ Nina Amenta∗ Dan A. Alcantara∗ Deboshmita Ghosh∗ Yong J. Kil∗ Eric Delson‡ Will Harcourt-Smith‡ F. James Rohlf§ Katherine St. John‡ Bernd Hamann ∗ ∗ Institute for Data Analysis and Visualization and Department of Computer Science, University of California, Davis ‡American Museum of Natural History and Lehman College of the City University of New York §Department of Ecology and Evolution, State University of New York, Stony Brook A BSTRACT 1 I NTRODUCTION We introduce a technique to visualize the gradual evolutionary Darwin’s theory of evolution was originally applied using morphol- change of the shapes of living things as a morph between known ogy - discrete qualitative features such as number of toes, and also three-dimensional shapes. Given geometric computer models of quantitative shape differences, such as elongation of a limb - to anatomical shapes for some collection of specimens - here the place species within the tree of life. While genomic sequence data skulls of the some of the extant members of a family of monkeys is now the basis of most phylogenies (evolutionary trees), morphol- - an evolutionary tree for the group implies a hypothesis about the ogy continues to be an essential part of evolutionary biology. One way in which the shape changed through time. We use a statisti- important reason is that the morphology of fossils, rather than com- cal model which expresses the value of some continuous variable at parisons between genomic data, provides our only direct evidence an internal point in the tree as a weighted average of the values at for extinct species. the leaves. The framework of geometric morphometrics can then be For instance, the shape of the skull is very important in the study used to define a shape-space, based on the correspondences of land- of human evolution, and that of our primate relations. Quantita- mark points on the surfaces, within which these weighted averages tive differences, such as shape of the brow ridge, are essential in can be realized as actual surfaces. defining the criteria for comparison between skulls. The idea that quantitative shape differences can be analyzed in terms of the trans- Our software provides tools for performing and visualizing such an formations required to “morph” one shape into another goes back analysis in three dimensions. Beginning with laser range scans of at least to d’Arcy Thompson’s classic 1917 book On Growth and crania, we use our landmark editor to interactively place landmark Form [24], from which Figure 1 is taken. points on the surface. We use these to compute a “tree-morph” that smoothly interpolates the shapes across the tree. Each intermediate shape in the morph is a linear combination of all of the input sur- faces. We create a surface model for an intermediate shape by warp- ing all the input meshes towards the correct shape and then blending them together. To do the blending, we compute a weighted average of their associated trivariate distance functions and then extract a surface from the resulting function. We implement this idea using the squared distance function, rather than the usual signed distance function, in a novel way. Figure 1: Spatial relationship between human, chimpanzee, and ba- boon skull, as envisioned by d’Arcy Thompson in 1917. The overall shape is recognizably similar, and the transformation between them CR Categories: I.4.10 [Image Processing and Computer Vision]: describes the shape difference. In modern morphometric studies, Image Representation—Morphological; I.3.8 [Computer Graph- the statistical analysis of the transformation is based on a matrix of ics]: Applications selected landmark points from the surfaces, while warps of the em- bedding space, which are illustrated here, are more often used for visualization of the results. Keywords: morphometrics, morphing, surface blending, merging, warping, distance fields, extremal surface The statistical analysis of geometric shape transformations is the ∗ e-mail:

,

, program of geometric morphometrics. In addition to evolutionary

,

,

biology, morphometric techniques are used widely in developmen- tal biology, medical image analysis, and other areas. Morphomet- rics defines “shape spaces” based on sets of homologous sets of landmark points on the input objects. The spaces in which statisti- cal analysis is generally done are linear spaces spanning the input shapes, so that the interpolating shapes in the space, such as those closed or oriented surfaces as input nor do we require any prepro- we use here, are linear combinations of the input shapes. cessing of the laser range scan surfaces. Using morphometrics, we have implemented a three-dimensional “tree-morph” visualizing the evolutionary changes implied by a given evolutionary tree. Surface meshes captured by a laser range 1.1 Geometric Morphometrics scanner from cranial specimens for extant (living) species appear at the leaves. The interior nodes and interior points of the edges of the evolutionary tree correspond to hypothetical ancestor species. We Geometric morphometrics is a branch of biostatistics dealing with realize these by computing synthetic surface meshes for the shapes the analysis of shape [4, 27, 1]. Scientists need to be able to de- at the internal nodes and at a dense set of points sampled along the fine and analyze statistically significant variables expressing bio- tree edges. We then visualize the morph by displaying the precom- logical shape. This task is difficult because the choice of what to puted meshes interactively, by sliding a cursor along tree edges, or measure and analyze affects the results. Rather than measuring spe- as an animation. cific distances, angles, and so on, the approach used in geometric mophometrics is to choose a discrete set of K homologous (or “cor- This specific project was chosen to drive the development of com- responding”) landmark points Li = {x1 , x2 , x3 , . . . , xK }, 1 ≤ i ≤ N puter tools for three-dimensional morphometric analysis and visu- on all N input object surfaces. The representation of the shape by alization. Our aim is to use visualization and interactive tools to its set of landmark points subsumes measurements of specific dis- support the kind of morphometric analysis that is currently done tances between landmarks, angles produced by three landmarks, in paleontology and other biological disciplines, producing images etc. A dense enough landmark point set provides an adequately which reflect the users’ existing theories about shape transforma- sampled representation of the shape. Statistical analysis of land- tion and which are therefore perceived as credible and relevant to mark point sets provides a method for making assertions such as the science. We do not propose in this initial project to change the “This fossil cranium resembles a macaque rather than a baboon,” current analytic techniques. But we do expect that as better tools more precise and quantifiable. make it possible to handle large amounts of three-dimensional data a greater emphasis on automation will be inevitable, and that the The sum-squared distance between two sets of landmark points is software framework we have developed here will make it easier to innovate. s K 2 We faced a number of design decisions and challenges in the i j D(L , L ) = ∑ Lin − Lnj (1) project. Morphometric techniques are based on user-defined sets n=1 of landmark points, which are assumed to be exactly homologous. Recent work in graphics has tended to instead to emphasize the au- Procrustes distance between them is the minimum of D(Li , L j ) tomatic determination of homology. This would be appropriate to over all rotations, scales, and translations of L j (with the scale of Li the extent that it reliably corresponds to biological homology, but normalized). Procrustes distance imposes a geometry on the space this is a terrific research question rather than an accepted technique of landmark configurations, forming a non-linear shape space. Un- to be applied. Our landmark editor supports existing practice by fa- fortunately, the pairwise alignments do not produce a multiple mu- cilitating the placement of large sets of landmarks, using automatic tual alignment of all the landmark sets; if we align La to Lb and Lc placement only very conservatively. This has had immediate and to Lb , we find that La is not aligned to Lc . The common practice is obvious impact. to compute an averaged consensus landmark point set and align all A second issue was producing the intermediate surfaces which lin- of the landmark point sets to this consensus configuration. Using early interpolate the input shapes, as defined in morphometric anal- a process called General Procrustes Alignment (GPA), a consen- yses. The multiple-alignment and interpolation procedures we use sus landmark set is chosen to minimize the total squared difference to handle the landmark point sets, while standard in morphometrics, between the aligned input landmark sets and the consensus config- differ from those common in computer graphics. They are designed uration [9]. Statistical analysis can then be performed on the data to minimize the error induced by forcing the shapes into a linear matrix formed by the aligned landmark sets. This is treated as a lin- space. We then needed to develop a procedure to produce linearly ear space, in which standard statistical techniques can be applied di- interpolated surfaces as well as linearly interpolated sets of land- rectly. This linear space is an approximation to the non-linear shape mark points. This is a different problem from standard morphing space defined by the Procrustes distances, and the choice of a good or blending. Besides having several inputs instead of just two (not consensus configuration is important to reduce the error caused by all methods generalize in an intuitive way), we wanted to preserve forcing the data into a linear space. properties of linear interpolation. In particular, we wanted the oper- ation to be commutative, meaning that, for instance, the surface one After GPA we can linearly interpolate our set of input configura- third of the way from A to B (2/3A + 1/3B) should be the same as tions, producing intermediate configurations which also lie in the the surface two-thirds of the way from B to A (1/3B + 2/3A); this approximating linear space. All of our intermediate landmark sets is not true for many existing two-input morphing methods. This in the tree-morph belong to such an interpolating linear space, with motivated our choice of a morphing method based on linear inter- the weights for the linear combinations determined as described in polation of implicit functions. Section 2.2. A third challenge was handling the somewhat messy scanned in- It is common practice in morphometrics to visualize an averaged put data, which does not even approximate a closed manifold sur- configuration by warping one of the input surfaces so that its land- face; there are large holes and only one side of the solid is captured marks coincide with the averaged configuration using a thin-plate in most areas, while in others both sides are very close together spline (TPS) warp. The TPS is favored because the resulting sur- or even self-intersecting. In the implict function interpolations, we face interpolates the landmarks (which are the “real” data being use the squared distance function to the surface rather than the more visualized), because it is the smoothest such warp, minimizing the usual signed distance function, from which we extract an extremal bending energy, which is related to the curvature, and because it is surface (defined in Section 2.5). This approach does not require fairly straightforward to compute. 1.2 Application to Primate Evolution The early steps follow the conventional course of a geometric mor- phometrics analysis, which is the “gold standard” for the scientists We use this well-accepted morphometric framework to define the in terms of how they model shape change. We relay on user-defined biologically “correct” interpolation of the skull shapes for five landmark points, which we accept as truly homologous. We focus species of Old World monkeys. The Old World monkeys evolved on making it easy for the user to define many landmarks, provid- in the same time and place as early humans, making them a partic- ing more data for the landmark-based statistical model. We use the ularly interesting group to study. There are many extinct species of GPA alignment procedure followed in morphometric analysis. Our Old World monkeys, known from fossils, so that there is a lot of in- warping step uses the TPS warp, which is well-accepted for the rea- teresting data related to their evolutionary history. Yet this history, sons discussed above. The final blending step meets two objectives. defined as the exact shape of the evolutionary tree, continues to be First, it carries through the principle of constructing the intermedi- a matter of controversy. ate shapes of the morph as a linear combination of the input shapes, by representing the shapes as trivariate functions for which linear We have used a laser range scanner to capture three-dimensional combinations are unambiguously well-defined. Second, our choice shapes of the crania of many species of Old World monkeys, both of trivariate function works well for the raw captured data meshes, extant and fossil, as part of a larger effort at the American Mu- which are not manifold and have holes, and which contain very thin seum of Natural History to develop a database of three-dimensional shell-like regions and occasional self-intersections which are diffi- primate morphology. Here, we use five of these surface models cult to morph. to develop and apply a method for visualizing morphometric esti- mates of skull-shape variation that is relevant to the evolution of the group. We selected five extant species sampling both subfami- 1.4 Related Work lies of Old World monkeys, and use a best-estimate for the evolu- tionary tree of the five species, derived from DNA sequence data, Existing geometric morphometric software has mainly focused on which is only available for the extant species. We then visualize the alignment and multivariate statistical analysis of specimens, what the sequence-based tree implies about the morphology of an- with less emphasis on either the landmark placement user inter- cient monkeys by interpolating the cranial shapes across the fixed face or on visualization. Morpheus [22], morphologika [18], and tree. Figure 2 shows the tree. the TPS suite of programs [19] are the packages most widely-used Visualizations of the intermediates (the hypothetical species at any by morphologists. interior point of the tree) are interesting in at least two ways. Sci- Placing landmark points on 3D specimens for morphometric anal- entists want to answer questions like, “Are the skull shapes pre- ysis is generally done using 3D contact digitizers on the actual ob- dicted by this model biologically plausible?” and “Where would jects where the collected points stored in a spreadsheet. This is this known fossil fit into the tree we hypothesize from genomic extremely tedious, so that landmark point sets consisting of tens data?” The visualization of the subset of the skull shape-space de- of points are typical. For virtual 3D images of specimens, such as fined by the tree helps to answer both kinds of questions. laser range surfaces or CT scans (ideal when the actual specimen cannot be obtained firsthand), there are generic software packages can be used, but these programs are not specialized for landmark 1.3 Technical Overview placement, so the process remains quite cumbersome. The interesting visualization problem of morphing primate skull Our goal is to produce a three-dimensional tree-morph visualiz- shapes across an evolutionary tree was first approached by Delson ing the evolutionary hypothesis presented by a specific evolutionary et al. [8] using the three-dimensional analog of the practice from tree using as input surface models captured by a laser range scan- morphometrics mentioned above, in which the transformation from ner and the transformations specified by the morphometric model. one shape to another determined by the landmark points is visual- The output of our procedure is a set of polygonal surface models, ized by warping one of the input surface models. This approach each one representing an intermediate shape corresponding to an has the drawback that the visualization of an intermediate produced interior point of the tree. Each of these intermediate models is a by warping one input surface is not the same as the visualization weighted average of the input models; they differ only in the choice produced by warping another instead. Our work improves on this of the weights. We developed a weighted average blending proce- approach in that we produce a single surface for each intermediate dure, applicable to rigid objects. It is not intended to handle inputs that represents the desired proportions of the input shapes. Also, in which the conformation varies as well as the shape (such as an since we can generate many more landmarks, we achieve a better arm bending at the elbow). Here is a brief description of each step; representation of the shape and its variation in the resulting model. more details are provided in the following sections. When the shapes are captured by computed tomography (CT) rather 1. Landmark specification: The user interactively places land- than laser range scans, the trivariate density functions for the differ- mark points at biologically meaningful locations providing ent specimens can be blended, after warping to align significant homologous points on each of the input specimens. features. This idea has been applied to visualizing the evolution of 2. Alignment and target computation: For each set of weights, toads by Hodges et al. [10]. The problem of merging similar sur- we produce a weighted average target configuration for the faces is replaced, in this case, with the problem of isosurfacing as landmark points, and aligne the input landmark sets to the the function is averaged across time, which is also non-trivial. target configuration, using GPA. Subsol et al. [23] produced an interesting morph between a CT scan 3. Warp: We compute a TPS warp taking the landmark points of a modern human skull and that of a fossil humaniod. Their goal on each input surface to the target configuration and use it to was to compare the two shapes using a deformation, and to demon- warp the entire input surface. strate some possible applications of three-dimensional geometry processing in paleontology. They establish correspondences using 4. Blend: We compute a distance function for each warped sur- automatically selected crest lines; the crest lines are noisy, the cor- face and take their weighted average. Extracting an extremal respondences between them are approximate, and they fail to cover surface from this function produces the output surface. many areas of the skull such as the brain case. So Subsol et al. Figure 2: The input surface meshes, from laser range scans of the crania of existing monkey species, are shown on the right-hand side at the leaves of this tree. The surface meshes at the internal nodes, placed at the estimated dates at which the species are hypothesized to have diverged, represent the skull shapes of the hypothetical ancestors as computed using our system. used them with a non-interpolating warp including a regularization landmark editor. The morphing method which formed their back component. This is a demonstration of the potential of possible end was more focused on efficiency and less on the accuracy of the automatic techniques, while we concentrate on the facilitation and intermediate shapes than some later approaches, including ours, and visualization of existing techniques which are actually used in sta- their method for color blending is somewhat similar to ours. Our tistical shape analysis. morphing method is most closely related to that of Levin, Cohen- Or and Solomovici [7] (and see also [25]). Their method uses the We also draw on methods known in computer graphics and visual- TPS to warp the surfaces so that they resemble each other closely. ization. We were inspired by one particularly relevant project [2], This nicely coincides with the common practice in morphometrics, in which a collection of full-body scans of humans was aligned to a where the TPS is favored. The surfaces are then merged by con- closed synthetic “base mesh.” The base mesh could then be warped verting each input surface into a signed distance function defined to resemble any of the inputs, or a linear combinations of the in- over a finite three-dimensional domain in the target space, taking a puts. This method produces a “space of human body shapes” useful weighted average of the functions, and extracting the zero-set of the in computer graphics, for instance for generating crowds of digital resulting function. This is appealing in our application since we can extras. In morphometrics, there is a strong emphasis on produc- take a weighted average of several input functions in a straightfor- ing results that are derived from the data rather than introduced for ward way. Their method works well for closed manifold input sur- computational convenience, so we wanted to avoid the synthetic faces, and reasonably well on our messy inputs except where there base mesh; also, we had no appropriate mesh to use. are self-intersecting surfaces or oppositely-oriented surfaces pass- ing very close to each other. We developed an alternative method Instead, we use a warp-and-merge paradigm to produce the inter- based on averaging the squared distance function, which produces mediate surface models. An early example of this paradigm is due a single-sheet output in such areas, with somewhat better results. to Lerios, Garfinkle and Levoy [12]. Their system included a user interface which allowed user to establish correspondences between Other morphing paradigms, which are better in the traditional curves and regions on the input models, similar to features of our graphics context where a surface A is morphed into another surface B, are not appropriate in our application. For example, the attrac- orientation of curves and the orientation and rotation of patches is tive level-set method of Breen and Whitaker [6] produces a morph shown with arrows, since it is easy to get these swapped (by placing by moving every point of surface A with velocities controlled by the control points in different orders on the two surfaces), and the the distance field of B. We do not see a straightforward way in user can re-orient a patch or curve to correct a mismatch without which this can be modified to produce a linear combination of sev- having to move the points. eral input surfaces. Also, the intermediate shape one-third of the way from A and B in the level-set morph is different from the in- Bookstein [5] introduced a method for optimizing the positions of termediate shape two-thirds of the way from B to A, so it cannot semi-landmarks on a curve or surface, minimizing the bending en- reasonably be interpreted as a linear interpolation. ergy of the induced TPS warps. We have implemented this method, and it seems to have minimal effect on our semi-landmarks, which are generally well-spaced to begin with. 2 S YSTEM D ETAILS We have also implemented a method that transfers landmark prim- itives from one surface to another semi-automatically. The user 2.1 Landmark Specification places at least four landmarks to produce a crude warp, which is then used to transfer the rest of the landmarks; the user then has to An essential part of the project was developing the interactive land- adjust their positions, but just transferring the overall configuration mark editor. A basic, but important, feature is that the homology be- greatly simplifies the experience and reduces errors. We can also tween the landmarks on different input surfaces is shown explicitly; export the landmark points, which allows existing morphometric with conventional methods, the user had to imply the homology by packages to use them. carefully placing landmarks in a specific order. In the landmark ed- itor, two surface meshes are shown at the same time in the main window. Figure 3 shows a screenshot. A dialog window shows the correspondences between pairs of landmark points; the windows are linked, so that when editing a correspondence in the dialog win- dow the selected points on the surface meshes are highlighted. Figure 4: A full set of 853 landmarks placed on one of the scanned crania. These were created using 45 are single points, 32 curves and 9 surface patches. Figure 3: Screen capture of our landmark editor. Two input meshes Using this interface, it is easy to create large sets of correctly cor- are shown in the large pane and upper-left pane while the two warped responding landmark and semi-landmark points (Figure 4). While models are overlayed in the lower-left pane. The yellow arrow indi- placing these was indeed tedious (it took our novice users about cates the selected patch orientation. three hours per skull), it would have been completely infeasible us- ing previous methods Points can of course be added, deleted and adjusted in any order. We show the surface normal as well as the point itself as the user The landmark editor is currently being used by the primatologists adjusts the landmark, which helps to place it exactly, especially on on our team for a separate research project investigating congruence high-curvature features. between joint surfaces in the primate skeleton, in which a grid of points are placed on the opposing joint surfaces. Landmarks have Not all shape differences can be captured by single points. To cap- been collected on laser range scan surfaces of over 80 primate lower ture curvature of an eye socket or the area of the brain case, we want limb-bone specimens, and results are being analyzed. The software to establish correspondences between curves and surface patches is greatly facilitating an otherwise lengthy and complex process. as well as points. Since morphometrics is based on the analysis of matrices of homologous points, in morphometrics this is done by distributing points, called semi-landmarks, across such features. We implement this by allowing the user to place the control points 2.2 Weights of B´ezier curves or patches on the surface. The system automati- cally generates a user-controlled number of semi-landmark points Each internal point of the tree corresponds to a set of consensus across the curve or patch, in a fixed order derived from the order of landmark points, which is a weighted average of the landmark the control points. The semi-landmarks are then projected onto the points at the leaves. The weights are determined using a princi- surface. The user establishes correspondences between the control ple known in evolutionary biology as squared-change parsimony: points two curves or two patches, implicitly defining the correspon- the integral of the squared change of a variable v (in this case, a sin- dences between all pairs of semi-landmarks on the primative. The gle landmark point coordinate) over the tree is minimized, within the constraints imposed by the values of v at the leaves. This princi- surfaces into this target space results in all surfaces being close to ple is sometimes used to estimate the structure of the evolutionary one another. tree from the values at the leaves [26]. Here instead we have the much easier problem in which the structure of the tree, including Conveniently, it is possible to get a closed-form solution for the the lengths of the branches, is fixed, and we just want to compute TPS warp, since it can be expressed as a weighted sum of radial the values of the variables at internal points in the tree. basis functions centered at the landmark points. The weights are determined by solving a linear system, which we do in a straight- A generalized least-square method for this problem was introduced forward manner. by Matrins and Hansen [14]. A covariance matrix for the values of v at the nodes is derived from the structure of the tree and used to weight a least-squares fit of the values at the internal nodes to the 2.5 Merging actual leaf data. The value v0 of v at the root is assumed to change randomly as it evolves through the tree. The value of v at some The final step is merging the N surfaces, which after the warp lie point in the tree is assumed to follow a Gaussian distribution, with very close together. This merge interpolates the surface details. mean v0 and and variance proportional to the distance from the root (this is the Brownian motion model of evolution). The covariance We initially implemented the method of [7], computing a signed for two points x and y in the tree is proportional to the distance from distance field around each surface, computing the weighted aver- the root to their least-common ancestor (lca); we assume the values age of those distance fields, and extracting the zero set. This works of v changed independently on the two separate paths from the lca well in most regions, but introduced holes in regions where both to x and y. sides of a thin sheet-like region of the object are captured in the input scans. If these thin sheets are not perfectly aligned by the This method is discussed in Rohlf [20]; his Equations 16 and 17 warp, holes appear in the output surface, even at arbitrarily high give the weights found by the least-squares fit as closed-form for- resolution (Figure 6). This occurs because even in the continuous mulas. We use them to assign, for a given internal node in the tree, case there are regions where the averaged function is positive over a set of weights w1 , . . . , wN for each of the N input surfaces. Along the whole neighborhood of the input surfaces and does not achieve each tree edge, we linearly interpolate the weights of the endpoints. zeros. These (and other artifacts near the boundaries) led us to con- sider an alternative approach, based on linearly interpolating the squared distance function rather than the signed distance function 2.3 Alignment and Target Configuration of the surface. Our new method is also somewhat less sensitive to discretization artifacts. Given N sets of homologous landmark points Li , 1 ≤ i ≤ N and a The squared distance function d i (x, y, z) associated with one input set of weights w1 , . . . , wN , the next step is Generalized Procrustes surface M i is zero at points of M i and positive elsewhere, and its Alignment (GPA) [9]. This is an iterative procedure that simulta- gradient at points of M i is the zero vector; M i lies at the bottom neously determines the positions of the landmarks on the output of the two-dimensional “valley” in the three-dimensional function surface and aligns the input landmark sets to this target configura- d i . The weighted average of the d i , d(x, y, z) = ∑ wi d i , is similar tion. to a squared distance function, in that there is still a distinct two- We begin by scaling each input set of landmark points so that the dimensional “valley” in the function. The points at the bottom of sum of the squared distances between allof the points and the center this valley have small but not zero values of d, and the gradient is of gravity is one. Thus, if gi = ∑K i exactly zero only at a set of discrete minimum points. We extract n=1 Ln /K is the center of gravity the surface tracing out the bottom of this valley. Notice that even at for point set i and the sum for each set is d i = ∑K i n=1 kLn − g k i 2 points where the averaged distance function fails to achieve a zero, i then the scaling factor for each point set relative to the first is s = the averaged squared distance function still has a valley. d i /d 1 , 1 ≤ i ≤ N. The scale variation can be re-introduced into the visualization at the end, as was done to make the images in Figure 2. We pick one of the input configurations of landmark points, arbi- trarily, as the first iteration of the target configuration. Then in each iteration, we align all of the input landmark sets to the proposed target configuration. We use Horn’s algorithm [11], which gives an efficient and closed form solution to the problem of finding the ro- tation and translation of an input landmark set minimizing the sum of squared distances between each landmark point and the homolo- gous point in the target configuration. After aligning all input land- mark sets, we compute the new target configuration by taking the Figure 5: Averaging multiple squared distance functions produces a weighted average of all of the homologous copies of each landmark function that is similar to a squared distance function, but generally point. We terminate the iterative process when the target configura- not zero anywhere. In one dimension, the squared distance function is a parabola in two-space, and the average of several is also a parabola. tion is “stable.” Fewer than ten iterations are typically needed. In higher dimensions the situation is similar, although complicated by the fact that the input surface normals do not match exactly. Taking a directional derivative in a direction v, roughly perpendicular to the 2.4 Warp desired surface, produces a signed function, with its zeros defining the bottom of the parabolic “valley.” The next step is to warp all surface models into the target configu- ration space using the thin-plate spline (TPS). The TPS warp is de- To extract a surface from a two-dimensional valley in a trivariate fined using an input set of landmarks and the target set of consensus function we use an extremal surface construction. Extremal sur- landmark points, and it brings the input landmarks into coincidence faces were introduced for extracting surfaces from noisy tensor with the target set. Taking the surface mesh through the same trans- fields [17], and they have been used recently for defining surfaces formation warps the surface into the target space. Warping all N from point clouds [3]. We consider the directional derivative of d in a direction v roughly perpendicular to the valley. Once we have a surface point p such that a sphere centered at x is tangent to the an appropriate vector field v, the directional derivative surface at p. The nearest foot-point is also the closest surface point. We modified the code to find not only the closest foot-point but ∂d ∇d · v also the second-closest foot-point as well. When the object forms g= = (2) ∂v 2 a thin sheet, or at a sharp corner, the second-closest foot-point is at nearly the same distance as the closest. We use the difference be- is a signed function, whose zero set is taken as the desired output tween the distance to the closest and the second-closest foot-points surface (Figure 5). Since the function g is locally nearly linear, as a measure of the difficulty of orienting v. For efficiency, we also extracting the output surface using marching cubes [13] works well modified the CPT code to only compute the squared distance field and produces smooth surfaces free of “jaggy” artifacts, except again around each input surface mesh M i to within a small distance α (we in regions where the input formed thin sheets (Figure 6). use 1/10 the largest dimension of the model). We use the CPT code to find the trivariate squared distance function d i for each M i and also the exact gradient ∇d i of the squared distance. Note that the gradient is the same as the unsigned distance function times two. A weighted average of all d i functions produces a single scalar trivari- ate function d, and, because of the linearity of the derivative opera- tor, averaging the ∇d i functions produces its exact gradient ∇d. Color is also averaged. Each triangle in the input mesh is assigned the average color of its three vertices (we could have interpolated the vertex colors across the triangle but the differences are negligi- ble at the grid resolution we are using). Every voxel closest to this triangle inherits this triangle color. The colors at the voxels are then averaged along with the distance functions. After the surface is ex- tracted, color is assigned to each surface vertex by interpolating the surrounding voxel colors. The resulting surfaces still have some extraneous parts. One issue is that g = 0 at ridges of d as well as at valleys, that is, near the the medial axis of the desired surface as well as the surface itself. We handle this simply by taking the single largest connected compo- nent of the output surface; this removes the medial axis components and also some small noise artifacts. We also delete any part of the Signed Unsigned surface where |g| is not nearly zero (zero crossings at transitions between large positive and large negative values of g). These occur Figure 6: Comparison of signed distance field blending, on the left, near boundaries. and our novel unsigned blending procedure on the right. In the first row, we use a high resolution voxel grid of 300 × 192 × 147. The Finally, at voxels where v is hard to define because the averaged signed distance field produces a merged surface with a hole behind gradients at the vertices point in very different directions, we do the eye socket, where two oppositely-oriented surface are close to not compute g directly. Instead, we compute g for the surrounding each other, even though none of the input surfaces had a hole. Our voxels and then interpolate the values, by averaging, to the empty method of averaging the squared distance functions produces the voxels. Note that this hole-filling would not be possible without a correct topology, although there is still some discretization noise in consistent orientation for v and hence a consistent sign for g. the difficult region. At lower resolution (in the second row, 150 × 96 × 73), our method is less sensitive to discretization artifacts. 3 R ESULTS To find a vector field v roughly perpendicular to the valley, we use the directions perpendicular to the input surfaces. Specifically, we average the unoriented gradient directions ∇d i of all of the input Figure 2 illustrates our results. The main point is that the synthetic squared distance fields d i ; by unoriented we mean the directions of skulls, created by averaging the input meshes, are virtually indistin- the lines supporting the vectors, so that the vectors v and −v are guishable from the original models. A video, including an anima- considered to be the same. Even in the difficult regions in which tion of the tree-morph and some examples of interaction with the the objects form thin sheets, the unoriented gradient directions are landmark editor, accompanies this paper. roughly consistent with each other (while the oriented gradient di- rections can point in opposite directions), so the averaged unori- The input surface meshes varied in size, from 797K to 433K tri- ented vector field is smooth near the surface. angles, except for the Papio model, for which only a 75K trian- gle mesh was available. Computing the trivariate distance function We then propagate a consistent orientation through v. Again, this from the input mesh is the most expensive part of the computation, has to be done carefully in the difficult regions in which the objects and this is roughly linear in the size of the input mesh. For the form thin sheets. We detect the difficult regions while extracting the animation, we simplified all of the meshes down to about 75K tri- distance function, as we describe in a moment. With this informa- angles, since more detail would not be resolved at video resolution. tion, we can propagate the orientation of v first in the easier regions Note that the trivariate distance function has to be computed for and then in the more difficult ones. each frame, since the warp for each frame is different. We used the full resolution input meshes for the figures. We compute the squared distance function di for each surface on a voxel grid using the closest point transform (CPT) code distributed We did our processing on four Intel 3.2GHZ Hyperthreaded work- by Mauch [16], which implements his robust computational- stations, each with 2GB of memory. The distance field for the high- geometry-based method [15]. This code finds for each voxel x the resolution meshes is computed in about 500 seconds per model on nearest foot-point on the polyhedral input surface; a foot-point is a voxel grid of size 300 × 192 × 147. Other processing, including the GPA, the TPS warp, and the extraction of the extremal surfaces, [3] Nina Amenta and Yong Joo Kil. Defining point-set surfaces. In Pro- required about an hour altogether and was minor compared to the ceedings of ACM SIGGRAPH, pages 264–270, 2004. time required to generate the distance functions. [4] F. L. Bookstein. Morphometric tools for landmark data: Geometry and Biology. Cambridge Univ. Press, New York, 1991. [5] F.L. Bookstein. Landmark methods for forms without landmarks: morphometrics of group differences in outline shape. Medical Image 4 D ISCUSSION AND FUTURE RESEARCH Analysis, 1(3):225–243, 1997. [6] D. Breen and R. Whitaker. A level-set approach for the metamorphasis This application raises a number of research questions that we are of solid models. IEEE Transactions on Visualization and Computer interested in pursuing. With respect to primate evolution, we plan Graphics, 7(2):173–192, 2001. [7] Daniel Cohen-Or, David Levin, and Amira Solomovici. Three- to compare the average ancestral shapes predicted by the statistical dimensional distance field metamorphosis. ACM Transactions on model and illustrated in this visualization with the shapes of known Graphics, 17:116–141, 1998. fossils, both visually and statistically. Integration of fossil evidence [8] E. Delson, D. Reddy, S. Frost, F.J. Rohlf, M. Friess, K. McNulty, with a trees such as ours, whose structure is inferred from genomic K. Baab, T. Capellini, and S. E. Hagell. 3d visualization of inferred sequence data on existing species, has to be based on morphological intermediates on a phylogenetic tree–applications in paleoanthropol- features. Visualizations such as these help paleontologists develop ogy (abstract). Amer. J. Phys. Anthropol. Suppl., 36:86–87, 2003. intuition about morphological change and encourage them to accept [9] J. C. Gower. Generalized procrustes analysis. Psychometrika, pages or reevaluate statistical models. 33–51, 1975. [10] W.L. Hodges, R. Reyes, Jr. T. Garland, and T. Rowe. Visualizing horn Generating landmarks automatically in a way that users would find evolution by morphing high-resolution x-ray ct images. In Sketches. sufficiently accurate and biologically meaningful is an important SIGGRAPH 2003, 2003. area for future research. As more data becomes available, the need [11] B.K.P. Horn. Closed-form solution of absolute orientation using unit for automation becomes more pressing. For instance, it would be quaternions. Journal of the Optical Society of America, 4(4):629–642, very helpful to be able to attract landmark points onto significant 1987. geometric features, especially ridges. More ambitiously, it would [12] A. Lerios, C.D. Garfinkle, and M. Levoy. Feature-based volume meta- be useful to be able to develop reliable surface correspondences us- morphosis. In Proceedings of ACM SIGGRAPH, pages 449–456, ing only a small number of landmarks, and hence transfer large sets 1995. of landmarks almost automatically. There has been some work on [13] William E. Lorensen and Harvey E. Cline. Marching cubes: A high this problem in the graphics community [21] and extending these resolution 3d surface construction algorithm. In Computer Graph- techniques to handle inputs that are not closed manifolds would be ics (Proceedings of SIGGRAPH 87), volume 21, pages 163–169, July helpful. 1987. [14] E.P. Martins and T.F. Hansen. Phylogenies and the comparative Finally now that we have a method for producing intermediate sur- method. American Naturalist, 149:646–667, 1997. faces, we can use it to construct a base mesh to use in a system [15] Sean Mauch. Efficient Algorithms for Solving Static Hamilton-Jacobi such as that of [2], which might produce more attractive or more Equations. PhD thesis, California Institute of Technology, Pasadena, efficient warps. California, 2003. [16] Sean Mauch. Closest point transform, 2004. http://www.acm.caltech.edu/ seanm/projects/cpt/cpt.html. [17] G. Medioni, M. S. Lee, and C. K. Tang. A Computational Framework ACKNOWLEDGEMENTS for Segmentation and Grouping. Elsevier Science B.V., 2000. [18] Paul O’Higgins and Nicholas Jones. Morphologika: Tools for shape We thank James Jones, a former U.C. Davis undergraduate, for analysis. writing the thin-plate spline code, and Tom Brunet, currently a [19] F. James Rohlf. tpsrelw: analysis of relative warps, 2004. [20] F.J. Rohlf. Comparative methods for the analysis of continuous vari- graduate student at the University of Wisconsin, for the code for ables: geometric interpretations. Evolution, 55(11):2143–2160, 2001. Horn’s algorithm. We thank Prof. Steve Frost for comments on the [21] J. Schreiner, A. Asirvatham, E. Praun, and H. Hoppe. Inter-surface manuscript, Lissa Tallman for beta-testing the landmark editor and mapping. In Proceedings of ACM SIGGRAPH, pages 870–877, 2004. John Liechty for sound recording. We thank Sean Mauch for mak- [22] Dennis E. Slice. Morpheus: Software for morphometric research. ing his CPT code available. We thank several anonymous referees [23] G. Subsol, B. Mafart, A. Silvestre, and M.A. de Lumley. 3d image for references and helpful suggestions. We thank the members of processing for the study of the evolution of the shape of the human the Visualization and Computer Graphics Research Group at the In- skull: presentation of the tools and preliminary results. In B. Ma- stitute for Data Analysis and Visualization (IDAV) at the University fart, H. Delingette, and G. Subsol, editors, Three-Dimensional Imag- of California, Davis. ing in Paleoanthropology and Prehistoric Archaeology, number 1049 in BAR International Series, pages 37–45. 2002. This work was supported by the National Science Foundation under [24] D’Arcy Wentworth Thompson. On Growth and Form. Dover Publi- contracts ACI 9624034 (CAREER Award) and CCR 0093378 (CA- cations, 1992 edition. REER Award), and the National Institutes of Health under contract [25] G. Turk and J. O’Brien. Shape transformation using variational im- P20 MH60975-06A2, funded by the National Institute of Mental plicit functions. In Proceedings of ACM SIGGRAPH, pages 335–342, Health and the National Science Foundation. 1999. [26] J. J. Wiens, editor. Phylogenetic analysis of morphological data. Smithsonian Institution Press, 2000. R EFERENCES [27] Miriam Zelditch, Donald Swiderski, David H Sheets, and William Fink. Geometric Morphometrics for Biologists. Academic Press, [1] D. C. Adams, F. J. Rohlf, and D. E. Slice. Geometric morphometrics: 2004. Ten years of progress following the ’revolution’. Italian Journal of Zoology, 71:5–16, 2004. [2] B. Allen, B. Curless, and Z. Popovic. The space of human body shapes: reconstruction and parameterization from range scans. In Pro- ceedings of ACM SIGGRAPH, pages 587–594, 2003.