Harmonic 1-form based skeleton extraction from examples

Graphical Models 71 (2009) 49–62

Contents lists available at ScienceDirect

Graphical Models journal homepage: www.elsevier.com/locate/gmod

Harmonic 1-form based skeleton extraction from examples Ying He *, Xian Xiao, Hock-Soon Seah School of Computer Engineering, Nanyang Technological University, Singapore

a r t i c l e

i n f o

Article history: Received 27 July 2008 Received in revised form 1 December 2008 Accepted 22 December 2008 Available online 24 January 2009

Keywords: Skeleton extraction Harmonic 1-form Deformation transfer

a b s t r a c t This paper presents a method to extract skeletons using examples. Our method is based on the observation that many deformations in real-world applications are isometric or near isometric. By taking advantage of the intrinsic property of harmonic 1-form, i.e., it is determined by the metric and independent of the resolution and embedding, our method can easily find a consistent mapping between the reference and example poses which can be in different resolutions and triangulations. We first construct the skeleton-like Reeb graph of a harmonic function defined on the given poses. Then by examining the changes of mean curvatures, we identify the initial locations of joints. Finally we refine the joint locations by solving a constrained optimization problem. We demonstrate the efficacy of the proposed framework by pose space deformation, skeleton transfer, shape segmentation and poseinvariant shape signature. Ó 2009 Elsevier Inc. All rights reserved.

1. Introduction Due to the increasing popularity of 3D laser scanners, highly detailed polygonal meshes are very common in industrial applications. Direct manipulating these models to produce faithful animation turns to be a tedious, timeconsuming and error-prone procedure. A common approach to deal with this problem is to use reduced models, such as skeletons [29], free-form lattices [46], surface markers or handlers [3], to ease the animation task. Among those reduced models, skeletons are perhaps the most intuitive and popular method in video games and computer animation industry. A well defined skeleton directly represents the anatomy structure of articulated characters, meanwhile it consists of highly reduced degrees of freedom compared with detailed models, and is straightforward for artists to manipulate, pose and animate characters. Extracting skeletons from 3D objects is challenging, especially in the context of finding exact centered skeletons. There exists large number of research efforts in find* Corresponding author. Fax: +65 67926559. E-mail addresses: [email protected] (Y. He), [email protected] (X. Xiao), [email protected] (H.-S. Seah). 1524-0703/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.gmod.2008.12.008

ing the skeletons of static models. Popular techniques include thinning [37,41], shape decomposition [24,30], Reeb graph [9,36], distance transformation [15], medial axis transform [2]. However, these methods are based on a single pose and rarely discuss skinning because some dynamic parameters, such as rigid bone transformations, joint locations, vertex weights, are usually difficult to obtain from a static pose. Therefore, it is desirable to use multiple poses instead of one. This paper focuses on the problem of example based skeleton extraction raised by Schaefer and Yuksel in [45]. We are strongly motivated by the observation: many real-world deformations are isometric or near isometric from the global point of view. For example, the horse has similar metric in the rest and running poses, i.e., the deformations are caused by bending instead of stretching. From the differential geometry point of view, if the two surfaces have similar metric but different embedding (appearance in R3 ), then their Gaussian curvatures are similar but mean curvatures are not (see Fig. 1). Usually, the deformations of articulated models, such as human and animals, occur at the joints. Therefore, the near-joint skins have similar Gaussian curvature but different mean curvature in the reference and example poses. On the other hand, the part of skins which are far from the joints can be considered

50

Y. He et al. / Graphical Models 71 (2009) 49–62

Fig. 2. Harmonic 1-form is an intrinsic property which is determined by the metric and is independent of the embedding and resolution. Although the rest and running pose have different triangulations and appearances, their metrics are similar, thus, the harmonic 1-forms are highly consistent. Harmonic 1-forms facilitate finding the one-to-one correspondence among different poses.

Fig. 1. Many deformations in real world applications are isometric or near isometric. As a result, the Gaussian curvatures do not change too much in various poses. However, mean curvature near the joints may change significantly. Please pay attention to the joints where the changes of mean curvatures are much larger than that of Gaussian curvatures.

as rigid, i.e., both the Gaussian curvature and mean curvature remain unchanged. Therefore, the curvatures can be used to identify and locate the joints. Although conceptually simple, the above idea is difficult to be directly applied to real-world examples. The reasons are twofold. First, the 3D shapes are usually acquired via 3D laser scanners. It is very common that the poses have quite different resolutions and triangulations. Finding the one-to-one correspondence between poses is non-trivial. Second, the commonly-used discrete operators to compute curvature are sensitive to the resolution and connectivity. The computation results of a given model in different resolutions may look quite different. Harmonic 1-form is a differential 1-form df defined on manifold M such that f : M ! R is harmonic. Harmonic 1form has many promising properties: – Harmonic 1-form is determined by the metric and is invariant under isometric transformation. – Harmonic 1-form is independent of surface representation. The surface in different resolution and triangulation has similar harmonic 1-form. – Harmonic 1-form is computationally efficient and robust as we only need to solve a linear system. – Symmetric harmonic 1-forms can be easily constructed on symmetric shapes. Fig. 2 shows that the harmonic 1-forms of the horse in two poses (with different triangulations) are highly consistent. Therefore, harmonic 1-form can serve as a

powerful tool for shape analysis, including the skeleton extraction. In this paper, we develop a robust algorithm for example based skeleton extraction. The key idea of the proposed approach is to first compute harmonic 1-forms of the reference and example poses. Due to the small changes of the metric, the isocurves of harmonic 1-forms across various poses are highly consistent. Then we use the isocurve based representation of the model and extract the skeleton-like Reeb graph of the harmonic functions. Next by examining the changes of mean curvatures, we identify the initial locations of joints. Finally we refine the joint locations by solving a constrained optimization problem. Compared to the existing work, the advantages of the proposed method include: – There is no restriction on the connectivity of the input meshes, i.e., the example poses could have very different triangulations and resolutions. – This method is automatic except that the user needs to choose a source point on the reference pose. – The number and locations of the joints and bones are totally determined by the dynamic geometry of the reference and example poses. – Since symmetric Dirichlet boundary values are set for symmetric shapes, the extracted skeletons also reflect the geometric symmetry of the input model. The remainder of this paper is organized as follows. Section 2 briefly reviews the previous work in skeletonization, Reeb graph, and harmonic 1-forms. Section 3 explains the algorithmic details for extracting skeletons from a set of example poses using harmonic 1-forms. Section 4 shows several applications of the proposed method, such as skeleton transfer and pose space deformation. Finally, we conclude the paper and briefly discuss the future research in Section 5.

Y. He et al. / Graphical Models 71 (2009) 49–62

2. Previous work 2.1. Skeleton driven deformation Animating characters driven by skeletons is widely applied though detailed description in academia is rare. Pioneer work can be found in [32] which introduced skeletal methods to animate hands and grasping behavior. Lewis et al. described standard method in industries thoroughly and obtained the name skeletal subspace deformation (SSD), which thereafter is also given by linear blending, soft skinning and single weight enveloping (SWE) [29]. Kavan and Zˇára identified that in SSD, linear blending of transformation matrix is ill defined for rotations and therefore proposes a novel method by blending associated quaternion instead [25]. Example based methods are proposed to correct underlying limitations of SSD. One type of methods is to correct collapsing parts by interpolating example shapes at run time [29,48,26,27,53]. The other category of approaches try to augment the SSD model by either introducing more parameters [55] or by adding extra joints to improve skeletal structure [35]. 2.2. Skeleton extraction Most of the existing work focus on extracting skeletons from a static pose of a model. Mortara and Patanè defined the affine-invariant skeletal representation for 3D shape matching [36]. Katz and Tal used fuzzy clustering method to decompose a shape and then extract the skeletons [24]. Lien et al. proposed an efficient and robust approach that simultaneously generates a hierarchical shape decomposition and a corresponding set of multi-resolution skeletons [30]. Theobalt et al. took a sequences of volume data to estimate skeleton through volume decomposition of motion data [50]. Oda et al. presented a framework of extracting skeleton interactively using geodesic distances [39]. Ma et al. constructed a RBF level set for 3D model and found the maximum of gradient descent for each point and then formed the skeleton [31]. Aujay et al. presented the harmonic skeleton in which anatomical information is used to enhance the skeletons [4]. They demonstrated that very realistic skeletons can be generated using harmonic functions as well as the prior of the anatomy of the articulated models. Dellas et al. presented an automatic method to extract the animation control skeleton of virtual humans, which relies on an a-priori knowledge of the human anatomy [11]. Baran and Popovic presented an automatic method to embed a generic skeleton to a wide range of articulated 3D characters [5]. Since it is usually difficult to estimate the dynamic behavior using the static pose alone, example-based techniques gain popularity in computer graphics in recent years [49,47,33,45]. In [22], James and Twigg demonstrated that the conventional skinning techniques can be extended to automatically skin deformable mesh animations. Instead of specifying the hierarchical kinematic skeleton, they estimated proxy bone transformations and vertex weights for deformable shape sequences. Recently, Schaefer and Yuksel proposed a novel method to extract hierarchical, rigid skeletons from example poses

51

[45]. They first defined ‘‘Rigid Error Functions” to find the best rigid transformation and used these error functions to estimate the transformations of the bones in the example poses. Then, they skinned the mesh by solving for vertex weights using a constrained optimization and bone influence maps. Finally, they determined the connectivity of the skeleton and the joint locations. Schaefer and Yuksel’s approach is capable to estimate the complete set of parameters for skeletal animation including bone transformation, skeletal hierarchy, joint location and vertex weights. However, their method requires that the reference and example pose have the same triangulation [45]. 2.3. Reeb graph Reeb graph is a powerful tool to analyze the topology of a shape [42]. It has wide applications in computer graphics, such as removing topological noise [54], skeleton extraction [9], shape abstraction and understanding [8,7]. Yoshihisa and Kunii proposed an algorithm to compute the Reeb graph in time Oðn2 Þ, where n is the number of edges of the mesh. Cole-McLaughlin et al. solved the problem in Oðn log nÞ time [10]. Recently, Pascucci et al. presented a robust online algorithm to compute Reeb graph [40] for extremely large dataset. 2.4. Harmonic 1-forms Given a scalar function f : M ! R defined on manifold M, harmonic 1-form df is a kind of differential 1-form such that f is harmonic, i.e., 4f ¼ 0. Harmonic 1-form plays a critical role in many applications of geometry processing. Gu and Yau pioneered global conformal parameterization using holomorphic 1-form, which can be decomposed to real and imaginary part and both parts are real harmonic 1-forms [17,23]. Holomorphic 1-form is also used to compute the affine structure of a given manifold, which is the key to construct manifold splines [16,19,20]. Guo et al. demonstrated that global surface conformal parameterization can be generalized to a more general setting of point based geometry by computing the holomorphic 1-form in a meshless manner [18]. Ni et al. used the harmonic Morse function to extract the topological structure of a surface [38]. Dong et al. [12] and Tong et al. [52] applied harmonic 1-forms in quadrilateral remeshing. 3. Skeletonization 3.1. Algorithm overview Our proposed method takes several poses of a given model as input. All poses are triangle meshes and could have very different triangulations and resolutions. Our algorithm runs in five successive stages, which are illustrated in Fig. 3 using the Armadillo model. (i) The user specifies a source point on the reference pose. We use the registration algorithm to map this source point to the other poses (Section 3.2). (ii) We compute the harmonic 1-form using the userspecified source point as constraints (Section 3.3).

52

Y. He et al. / Graphical Models 71 (2009) 49–62

Fig. 3. Algorithm pipeline. (a) The user selects a source point on the reference pose. Usually, this point locates on the head with salient features. (b) We compute the harmonic 1-form on all poses. Note that the given rest and example poses have the same triangulation. In order to demonstrate our technique, we intentionally remesh each pose in different resolution. Two of eleven poses are shown here. Since the metric are similar in various poses, the harmonic 1-form are highly consistent. Since the Dirichlet boundary constraints are symmetric on the armadillo model, the harmonic 1-forms, as expected, also reflect the symmetric of the shape. (c) Then, we construct the skeleton like Reeb graph of the harmonic function. (d) Next, we identify the joints by finding the isocurves where the mean curvature changes significantly compared to the reference and other poses. The color on the isocurves illustrates the average of mean curvature. (e) Finally, we find the optimal location of the joints and bone lengths by solving a constrained optimization problem. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(iii) We construct the skeleton-like Reeb graph of the harmonic function (Section 3.4). (iv) We extract the isocurves of the harmonic function and compute the average mean curvature for each isocurve. By comparing the mean curvature for all poses, we can identify the joints which are closely related to the isocurves whose mean curvatures are changed significantly (Section 3.5). (v) We optimize the joint locations and bone lengths by minimizing the distance between the skeletons and the Reeb graphs (Section 3.6).

3.2. Setting the source point Similar to [4], a source point is needed to set the boundary value for the harmonic function. This user-specified source point is one of the end points of the skeleton-like Reeb graph (see Fig. 3c). Thus, it is desirable to locate the source point such that there is a joint nearby. Furthermore, in our algorithm, the user only needs to specify the source point on one pose, then we use the registration algorithm to map the source point and its local neighborhood to other poses. To avoid the ambiguity of the registration, we require that the source point locates at the unique salient feature of the given model [13], e.g., the nose of the Horse and Armadillo model (see Figs. 3a and 5a). Note that if the input model is symmetric, the source point must locate on the symmetric plane to induce the symmetric harmonic 1forms. 3.3. Computing harmonic 1-forms Given a surface M with user-specified constraint vertices, vi , i ¼ 1; . . . ; k, define a scalar function f : M ! R, such

Fig. 4. Harmonic 1-form on genus two surface. The surface is cut along the red curve c in (a). Two boundaries are denoted by cþ and c . The Dirichlet boundary value of cþ and c are set to be +1 and 0, respectively. Note that symmetric harmonic 1-forms can be easily constructed on symmetric shape if the boundary constraints reflect the symmetric of the input model. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

that f ðvi Þ ¼ ai , i ¼ 1; . . . ; k, where ai are the user-specified constraints. The function f is harmonic if it minimizes the harmonic energy,

Eðf Þ ¼

Z

j 5 f j2

ð1Þ

M

Then 5f , the gradient of f, is a harmonic 1-form on M. The harmonic 1-form, denoted by x, is defined for each edge ½v; w, i.e., xð½v; wÞ ¼ df ð½v; wÞ ¼ f ðwÞ  f ðvÞ, which satisfies the following equation [17]:

4f ðvÞ ¼

X

kv;w xð½v; wÞ ¼ 0

ð2Þ

½v;w2M

The commonly-used weights kv;w include uniform weights and cotan weights. Harmonic 1-form can be computed either by solving the Laplace equation directly or finding the minimizer of the harmonic energy. Fig. 4 illustrates the harmonic 1-form on genus two surface. Specifying the boundary value is critical to the extrema of the harmonic function. Various methods are used in the

Y. He et al. / Graphical Models 71 (2009) 49–62

53

Fig. 5. Computing the harmonic 1-forms. The global minimal point vsource is drawn in red in (a). By solving the Poisson equation, 4pðvÞ ¼ HðvÞ and pðvsource Þ ¼ 0, we obtain a smooth function p which mimics the curvature of the surface shown in (b). Then we identify the local maximal points (four on feet and one on the tail tip) of p, which are shown in green in (c). Finally, the values of the extrema of the Poisson scalar field p are used as the boundary constraints of the Laplace equation. The resulting harmonic 1-form is shown in (d). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

previous work. Tierny et al. chose two points on the mesh which are the most distant [51]. Aujay et al. tried to find vertices v such that the distance dðvsource ; vÞ on the mesh is locally maximum [4]. Dong et al. suggested to solve a Poisson equation which mimics the curvature of the input mesh and then find the local extrema as the boundary vertices of the harmonic equation [12]. In this paper, we follow Dong et al.’s method [12] to automatically determine the locations and values of the boundary constraints once the source point vsource is given. For each pose, there is a unique global minimal source point vsource (see Fig. 5a) which is computed in Section 3.2. We solve the following Poisson equation as suggested in [12]:

4pðvÞ ¼ HðvÞ and pðvsource Þ ¼ 0

ð3Þ

where HðvÞ is the mean curvature of vertex v. The goal of solving the above Poisson equation is to find a scalar field p which mimics the curvature [12]. After clustering the local extrema points which are very close to each other, the representatives of each cluster are identified as the boundary points and the function values are used as the boundary values. In order to make the harmonic 1-form of the reference and example poses as consistent as possible, we apply the same Dirichlet boundary constraints across all poses. Note that many articulated models are symmetric so that to keep the symmetry of harmonic 1-forms, we also require the boundary values of symmetric points to be identical. For example, we compute the geodesic between

the source point and the local extreme points. If two or more extreme points have similar distances, we apply the same boundary value to them. This scheme guarantees the resulting harmonic 1-forms reflect the symmetric of the input shapes. Fig. 5 illustrates the harmonic 1-form on one of the example pose. 3.4. Constructing the skeleton-like reeb graph Reeb graph is a powerful tool when dealing with topological skeletons [42]. For a real-valued smooth function f : M ! R on M, the points whose derivatives of f vanishes are called critical points of f. The Reeb graph of f is a graph whose nodes corresponds to these critical points and encodes the connectivity between them. Note that the leaves of Reeb graph exactly match the local maximal and minimal of f. Building the graph. Reeb graph can be constructed by contracting the connected components of the isocurves (level sets) of f to a point. Let n be the user-specified number of samples of the isovalues (200 in our experiments). We first normalize the function value of f to the unit interval ½0; 1 and then uniformly sample isovalues, fi ¼ ni , i ¼ 1; . . . ; n. Next we scan every critical point p of f, rf ðpÞ ¼ 0. If f ðpÞ is not sampled, we insert f ðpÞ into the isovalue set. Denote ck ðfi Þ the isocurve with isovalue fi on kth pose. k ¼ 0 refers to the reference pose. (Note that each isovalue fi may have more than one isocurve. We will explain eliminating the ambiguity later.) For each isocurve,

54

Y. He et al. / Graphical Models 71 (2009) 49–62

we compute its center as the representative. Then we use the sweep algorithm to connect these representatives and form the Reeb graph [10]. Topological filter. Once the Reeb graphs of reference and example poses are constructed, we first identify the key nodes (valence – 2). Due to numerical error and/or metric changes, multiple key nodes may occur in locations where there should only be a single key node. Thus, we cluster key nodes such that the Reeb graphs of the reference and example poses have the same topological structure, i.e., the same number of key nodes (see Fig. 6). Eliminating the ambiguity. Note that due to the symmetry of many articulated models, such as human and animals, two or more isocurves with similar geometry may have the same isovalue, e.g., the Armadillo’s elbows. Thus, we then need to distinguish these isocurves with the same isovalues. In this paper, we take advantages of the articulated models which usually have some key points in the Reeb graph. Since the user picks the source point on the head, the first key point connecting to the source is the neck. So we can identify the spine easily. To further classify a branch into left and right arms (legs), we compute the sum of the signed distance from the points on the branch to the symmetry plane passing the spine. If the sum is positive, the branch is identified as a left arm or leg depending on the distance to the source point. Otherwise, the branch is identified as the right arm or leg. Note that the sign determining left and right is not important. The user can freely change the sign to swap left arm/leg and right arm/leg. Finding the correspondence among poses. Next we want to find the one-to-one correspondence between reference pose and example pose. Note that we don’t need to find such a mapping between the vertices of two poses which are usually in different triangulation. Instead we will find the mapping for isocurves. Isocurves are the level sets of the harmonic function f, which are independent of the resolution and embedding. Therefore, even though the reference and example poses have very different resolutions, their isocurves are highly consistent if the metrics are similar. However, in practice, there always exist some cases, especially the poses created by users instead of acquired by 3D scanners, where the deformations causes large changes of the metrics (although they should NOT change

too much). As a result, the harmonic 1-forms and the corresponding isocurves may differ a little bit. In order to find the best one-to-one correspondence of isocurves between the reference and example pose, we first compute the integral of Gaussian curvature kG along each isocurve i.e., R k dl. Then for every isocurve cik in the example pose ck ðfi Þ G k, k P 1, we find the isocurve in the reference pose with the same isovalue c0 ðfi Þ. If both isocurves have similar integral of Gaussian curvatures, we find the ideal mapping. Otherwise, we perturb fi of the example pose in a small neighborhood and find the isovalue/isocurve such that the difference of the integral of Gaussian curvatures is minimal. As shown in Fig. 8, the Reeb graphs after postprocessing are highly consistent. 3.5. Identifying the joints As explained above, many real-world deformations of articulated models are isometry or near isometry. It is well known that Gaussian curvature is completely determined by the metric, therefore, the Gaussian curvature of the reference and example poses remain almost unchanged. Mean curvature, however, is related to the embedding, and could be used to identify the joints. We observe that compared to the reference poses, the embedding of the joint-related skins could vary significantly in example poses. Thus, we use the average mean curvature

R  (R  ci kH dl k dl ci0 H   Rk R DðiÞ ¼ max     ci dl dl  ci k

0

 )   k ¼ 1; 2; . . . 

ð4Þ

to measure the difference of the embedding of isocurve ci between the reference pose and example poses. If DðiÞ reaches its local maximal, we mark the center of isocurve ci as a joint candidate for every pose. We should also point out that for the symmetric shapes, the harmonic 1-forms are also symmetric due to the symmetric Dirichlet boundary values are specified in Section 3.3. Therefore, for symmetric parts, e.g., the Armadillo’s arms, if one joint is identified, then the corresponding isocurves on the symmetric parts are marked as joints too. However, in our experiments, we experienced that using the above local maximal condition alone, some extra joints may be identified. For example, due to the muscle bulge in the example pose, the metrics in the bulged areas

Fig. 6. Topological filter on the Reeb graph. (a) One key node (valence – 2, drawn in red) exists in the reference pose. (b) Due to the small changes of metric in the example pose, two key nodes occur in nearby locations where there should only be a single key node. (c) We remove one key node by changing the local connectivity of the Reeb graph. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Y. He et al. / Graphical Models 71 (2009) 49–62

55

Fig. 7. Identifying the joints. Four of 16 poses are shown here. Row 1: The color on each isocurve illustrates the changes of the mean curvature between the reference and running poses. Row 2: By finding the local maximal of mean curvature changes, we can identify the joint candidates (green points) on the Reeb Graph. Note that the joints are closely related to the isocurves whose mean curvatures change significantly. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. Constructing the skeleton-like Reeb graph. (a) Isocurves of harmonic function f on the surfaces. (b) The red points are the critical points of the harmonic function f. The Reeb graphs of f on different poses are very consistent. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

change a little bit larger than the rigid areas (see Fig. 10). As a result, an extra joint is identified. To remove these unnecessary joints, we use the following angle constraint. For any valence-2 joint, we compute the smallest angle between the two incident bones in ALL poses. If the angle is greater than the user-specified threshold (160° in our implementation), this joint is considered as extra joint. As expected, this joint detection approach is highly related to the examples. The more poses involved, the more accuracy and robustness of the joints to be identified. Sixteen poses of horse model are used in our experiment. Fig. 9 shows the maximal mean curvature difference max DðiÞ of head and leg. Obviously, the joints are closely related to the isocurves whose mean curvature differs sig-

nificantly to the reference pose. Fig. 7 shows the detected joints on example poses. 3.6. Optimizing the joint locations Section 3.5 identifies the number of joints in the reference and example poses. However, the joints just locate on the Reeb graph. Although encoding the topological structures nicely, the skeleton-like Reeb graphs are not the anatomical skeletons. Therefore, we need to estimate the joint locations using examples. Assume m joints are identified using the above method. Denote Jjk the jth joint, and Blk the length of lth bone in kth pose. We enforce that the bone lengths do not change in the deformations. For an isocurve cik in the kth pose, denote

56

Y. He et al. / Graphical Models 71 (2009) 49–62

dðcik ; J 1k ; . . . ; J m k Þ the sum of distance between the point on the isocurve and the skeleton defined by the joints fJ jk gm j¼1 . Then, we also require that such distances are also invariant in all poses too. These constraints lead to the following constrained optimization problem:

arg min Jj k

2 X X  1 m 1 m  dðcik ; Jk ; . . . ; J k Þ  dðci0 ; J 0 ; . . . ; J 0 Þ kP1

ð5Þ

i

subject to the constraints

Bl0 ¼ Blk

Fig. 9. The column graphs show the maximal mean curvature difference DðiÞ of the head and one leg. The red bars are the local maximums which correspond to joints (green points). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

8l; k

ð6Þ

The objective function ensures that the relationship between the skin and bone is as rigid as possible. The constraints ensure that the bone length remains unchanged in all poses. Figs. 11 and 12 show the optimized skeletons of the horse and lion models, respectively. We tested our program on a workstation with 2.6 GHz CPU and 3 G memory. The statistics of the test examples are shown in Table 1. Note that no source point is specified for the Arm model. Instead, the Dirichlet boundary conditions are set as 1 and 0 on the two boundaries. 4. Applications 4.1. Cross parameterization and skeleton transfer Skeleton transfer is the direct applications of the proposed method. In this section, we demonstrate transferring the dynamic skeletons from the horse to the cat (see Fig. 13). We pick a source point on the cat’s nose. Then we compute the cross parameterization using harmonic 1-forms. To make the harmonic 1-forms of horse and cat as consistent as possible, we use the same boundary value conditions, i.e., the harmonic function values of the cat’s nose, feet and tail are exactly the same as the corresponding values of the horse model. Once the harmonic 1-forms are ready, we can easily build the oneto-one correspondence between the isocurves of horse and cat. Since each joint is associated with a unique isocurve, the skeleton of horse can be easily transferred to the cat. 4.2. Pose space deformation

Fig. 10. The muscle bulge may cause large change of metrics on the example poses (see Rows 1 and 2). As a result, some extra joints will be identified (see Row 3). To avoid these extra joints, we use the angle constraint to test each joint. If the angles between two incident bones are greater than the user-specified threshold in ALL poses, this joint is identified as extra joint and then discarded (see Row 4).

Pose space deformation (PSD) provides one solution for problems induced by common underlying skinning methods such as SSD. It ‘‘corrects” collapsing parts of deformed models on the fly by interpolating well sculpted examples. Most literatures on PSD assume a single example as a pair of geometric mesh with associated skeleton [29,48,26,27]. Building a skeleton and rigging it to polygonal meshes are usually tedious and time-consuming to artists, not to mention posing it to align deformed models with given examples. In previous sections, the extracted skeletons only compose of joint positions in world space. To employ the skeletons in real productions, it is desirable to compute joint transformations as parameters of pose space. Although an example mesh can be expressed as a reference mesh with a group of affine transformations [22], animators prefer

Y. He et al. / Graphical Models 71 (2009) 49–62

57

Fig. 11. Optimizing the joint locations by solving a constrained optimization problem.

Fig. 12. Extracted skeleton of the lion model (three of eight poses are shown).

Table 1 Statistics of test examples. N p , number of poses; Nt , average number of triangles in each pose (note that we intentionally remesh the given meshes to get the different triangulation of each pose); T s , time for mapping the source point from the reference pose to example poses; T h , time for computing harmonic 1-forms; T r , time for constructing the Reeb graph; T j , time for identifying the joints; T o , time for optimizing the joint locations. Time are measured in seconds. Object

Np

Nt

Ts

Th

Tr

Tj

To

Arm Horse Armadillo Lion

4 16 6 8

3k 20k 60k 12k

NA 10 4 3

12 119 270 50

6 51 23 20

1 10 8 4

3 40 30 23

more meaningful parameters to control such as joint rotations or translations. In this paper, instead of computing pure rotations [28] or rotations and translations [34], we parameterize example pose space as rotations and scaling

values based on the observation that scaling joints can capture more accurate deformations than translations. This strategy can also be justified by the fact that articulate deformation is more likely resulted from muscle flexion while bone length should be fixed. We first align the source joint of reference mesh to one example mesh. And then starting from root joint, rotation angles and scaling values for each joint will be recursively computed in a manner of width first search. Given joint positions of a pair of parent and child joint in reference pose and one example pose, two bones in respective reference and example pose form a quaternion that rotates reference bone into example bone. Converting this quaternion to euler rotations in the local coordinate system of one parent joint is trivial. The offset between rotated child joint in reference pose and the child in example pose can serve as the amount of joint scaling. Once pose space parameters of each example mesh are ob-

Fig. 13. Cross parameterization and skeleton transfer using harmonic 1-form.

58

Y. He et al. / Graphical Models 71 (2009) 49–62

tained, we can perform PSD with the method described in [55]. Fig. 14 shows that additional poses can be easily created using PSD. 4.3. Skeleton-driven shape segmentation Shape segmentation aims to decompose 3D surfaces into its ‘‘meaningful” components that a human being can intuitively perceives as distinct and logical parts of the shape. Among various approaches, skeletons have proven to be a very effective tool in shape segmentation since skeleton itself is a compact and meaningful shape descriptor [43]. In our current framework, we extract the skeletons from a large number of given poses. With the aid of extracted skeleton, we can segment the various poses in a consistent fashion. Our segmentation approach is inspired by the skin attachment method [5] in which the vertex weights are computed by solving a heat equilibrium equation. In order to get a meaningful segmentation, we need to know how much a bone could control the attached skins. Given a mesh M and extracted skeleton, let Bi denote the ith bone. Similar to [5], we solve the following PDE for every bone Bi ,

4wi þ Hðpi  wi Þ ¼ 0 where wi contains the weights of each vertex in M to bone Bi , pi is a vector with pij ¼ 1 if the nearest bone to vertex j is i and 0 otherwise, H is a diagonal matrix which measures the shortest distance from a vertex v to the nearest bone. The computed weights wi can be used in the linear blend skinning (LBS) method. Intuitively speaking, the weights wi measure how much each bone transform affects each vertex. Note that the computed weight wi is a smooth scalar function defined on M. We normalize the weights to range ½0; 1 and then trace the isocurve with the userspecified isovalue (0:5 in our experiments). The extracted isocurves serve as the cut locus for shape segmentation. Fig. 15 illustrates the proposed shape segmentation approach. Note that our extracted skeleton is the minimizer of optimization Eq. (5) where the shortest distances from a vertex to the nearest bone do not change across various poses, H and pi are also invariant to the poses. Therefore, the computed weights wi is independent of the poses. As a result, the proposed shape segmentation is poseinvariant.

Fig. 14. New poses construction using the extracted skeletons.

Fig. 15. Skeleton-driven shape segmentation. Given a mesh (a) and extracted skeleton (b), we first compute the weights wi for each bone. The weight wi is a smooth scalar functions defined on M (see (c)), the closer of a vertex to the bone, the larger the value of the weight. Then we trace the isocurve of the scalar function with the user-specified isovalue (0.5 in our experiments) as shown in (d). The segmentation is obtained by merging the regions of the userspecified joints (usually the degree-2 joints) (see (e)). Note that the weight function wi is pose-invariant, so is the segmentation (see (f) and (g)).

Y. He et al. / Graphical Models 71 (2009) 49–62

59

4.4. Pose-invariant shape signature

M/ ¼ K T  K orig

Shape signature plays an important role in shape analysis and retrieval. A shape signature captures the essential features of the shape and does not fully represent the shape. Thus, it is impossible to reconstruct the shape back from it. The articulated models, such as human or animal, may come in many different poses. Although a pose may look significantly different than the other, they represent the same object. The shape signature for articulated model is required to tolerate rotation, scaling, translation and deformation. Therefore, features depending on orientation, size and location can not be employed as shape signatures. Some examples of pose-invariant shape signatures include local-diameter function/centricity function [14], GPS embedding[44], spectral embedding [21], conformal factor[6]. In this subsection, we develop a new pose-invariant shape signature using isocurves of the harmonic function. In [6], Ben-Chen and Gotsman computed the conformal factor / by solving the following Poisson equation:

where K orig is the discrete Gaussian curvature of M in the original metric and K T is the constant target Gaussian curvature. The original and target Gaussian curvatures satisfy P P T orig ¼ 2pv, the Gauss–Bonnet Theorem, v2M K ¼ v2M K where v is the Euler characteristic. Note that K T , K orig and Laplace–Beltrami operator in the above equation are intrinsic, so is the computed conformal factor /. Therefore, it can be used for pose-invariant signature for articulated models. In [6], Ben-Chen and Gotsman used the histogram of the conformal factor / as the shape signature. However, conformal factor histogram alone only characterizes the distribution of the value of /, many important spatial information are lost. To fully utilize the computed /, we propose two new signatures, Normalized Length (NL) and Integration of Gaussian Curvature (IGC). Given a scalar function / : M ! R (/ can be either the harmonic function computed in Section 3.3 or the conformal factor computed by [6]), we first normalize the function value of / to the unit interval ½0; 1 and then

Fig. 16. Normalized Length (NL) and Integration of Gaussian Curvature (IGC) compared to Conformal Factor (CF) [6]. (a) shows the harmonic 1-form. (b) and (c) show the histograms of NL and IGC of the harmonic 1-form respectively. (d) shows the histogram of conformal factor.NL and IGC are more expressive regarding the distinction between shapes and more robust to pose changes.

60

Y. He et al. / Graphical Models 71 (2009) 49–62

gration of Gaussian curvature (IGC) along the isocurves in C i , i.e.,

IGCð/i Þ ¼

X

R

j

ci

R

j

ci 2C i

Fig. 17. Precision/recall graphs of Conformal Factor (CF [6]) and the proposed NL/IGC shape signatures for the SHREC benchmark.

uniformly sample isovalues, /i ¼ ni , i ¼ 1; . . . ; n, as shown in Row 1 of Fig. 7. Note that two or more isocurves may correspond to the same isovalues. Let C i ¼ fc1i ; . . . ; cki g denote the set of isocurves whose isovalue equals /i . Then we construct two functions for each isovalue, /i . The first function measures the normalized length (NL) of the isocurves in C i , i.e.,

NLð/i Þ ¼

1X j lðci Þ d j

ð7Þ

ci 2C i

where lðcji Þ is the length of curve cji , d is the main diagonal of the shape M. The second function characterizes the inte-

kG dl j

ci

dl

ð8Þ

Then our shape signature is defined as two histograms of NL and IGC. As shown in Fig. 16, IGC and NL are more expressive regarding the distinction between shapes and more robust to pose changes. Although the model is represented in many different poses, the shape signatures are highly consistent. We also apply the shape signature to different models with similar poses. As shown in Fig. 18, models of the same class have very similar histogram, and different models may have very distinctive histogram. To show the effectiveness of the proposed shape signature, we use the standard precision/recall graph for the SHREC benchmark [1]. Fig. 17 shows the precision/recall graph for our retrieval experiments. 5. Conclusions In this paper, we presented an approach to example based skeleton extraction. The proposed method is automatic except that the user needs to specify a source point on the reference pose. Due to many promising properties of harmonic 1-forms, such as resolution and embedding independence, we constructed a highly consistent one-toone correspondence between reference and example

Fig. 18. Shape signature of different models in various poses. It can be seen that models of the same class have very similar histogram, and different models may have very distinctive histogram.

Y. He et al. / Graphical Models 71 (2009) 49–62

poses. Therefore, even though the input models are in different resolutions, the extracted skeletons are highly consistent. The number and locations of the joints are totally determined by the dynamic geometry of the reference and example poses. Thus, the more poses are involved, the more accurate and robust skeletons we will obtain. Furthermore, for symmetric shapes, the extracted skeletons can also reflect the geometric symmetry of the input model. We demonstrated the proposed framework is effective for several applications such as skeleton transfer, pose space deformation, shape segmentation and pose-invariant shape signature. The limitation of the proposed method exists and demand further improvement in the future. Our current method only applies to the model which we can eliminate the ambiguity using the distance from critical points to the source point (usually on the head). However, this method does not work for octopus which has equal sized arms. In the near future, we plan to improve our method to apply to a wide range of objects. Acknowledgements This work was supported in part by the NTU SUG19/06, AcRF RG69/07, NRF2008IDM-IDM-004-006. The models are courtesy of Robert Sumner and Jovan Popovic (Horse) and Shin Yoshizawa (Armadillo). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.gmod.2008. 12.008. References [1] . [2] N. Amenta, S. Choi, R.K. Kolluri, The power crust, in: Proceedings of the Sixth ACM Symposium on Solid modeling and Applications, 2001, pp. 249–266. [3] O.K.-C. Au, H. Fu, C.-L. Tai, D. Cohen-Or, Handle-aware isolines for scalable shape editing, ACM Trans. Graph. 26 (3) (2007) 83. [4] G. Aujay, F. Hétroy, F. Lazarus, C. Depraz, Harmonic skeleton for realistic character animation, in: Symposium on Computer Animation, 2007, pp. 151–160. [5] I. Baran, J. Popovic, Automatic rigging and animation of 3d characters, ACM Trans. Graph. 26 (3) (2007) 72. [6] M. Ben-Chen, C. Gotsman, Characterizing shape using conformal factors, in: Proceedings of Eurographics Workshop on Shape Retrieval, 2008. [7] S. Biasotti, B. Falcidieno, M. Spagnuolo, Extended reeb graphs for surface understanding and description, in: DGCI, 2000, pp. 185– 197. [8] S. Biasotti, B. Falcidieno, M. Spagnuolo, Shape abstraction using computational topology techniques, in: From Geometric Modeling to shape Modeling, 2002, pp. 209–222. [9] S. Biasotti, S. Marini, M. Mortara, G. Patanè, An overview on properties and efficacy of topological skeletons in shape modelling, in: Shape Modeling International, vol. 297, 2003, pp. 245–256. [10] K. Cole-McLaughlin, H. Edelsbrunner, J. Harer, V. Natarajan, V. Pascucci, Loops in reeb graphs of 2-manifolds, in: Proceedings of the Nineteenth Annual Symposium on Computational Geometry, 2003. pp. 344–350. [11] F. Dellas, L. Moccozet, N. Magnenat-Thalmann, M. Mortara, G. Patanè, M. Spagnuolo, B. Falcidieno, Knowledge-based extraction of control skeletons for animation, in: Shape Modeling International, 2007, pp. 51–60.

61

[12] S. Dong, S. Kircher, M. Garland, Harmonic functions for quadrilateral remeshing of arbitrary manifolds, Comput. Aided Geom. Des. 22 (5) (2005) 392–423. [13] R. Gal, D. Cohen-Or, Salient geometric features for partial shape matching and similarity, ACM Trans. Graph. 25 (1) (2006) 130–150. [14] R. Gal, A. Shamir, D. Cohen-Or, Pose-oblivious shape signature, IEEE Trans. Vis. Comput. Graph. 13 (2) (2007) 261–271. [15] Y. Ge, J.M. Fitzpatrick, On the generation of skeletons from discrete euclidean distance maps, IEEE Trans. Pattern Anal. Mach. Intell. 18 (11) (1996) 1055–1066. [16] X. Gu, Y. He, H. Qin, Manifold splines, in: ACM Symposium on Solid and Physical Modeling, 2005, pp. 27–38. [17] X. Gu, S.-T. Yau, Global conformal surface parameterization. in: Proceedings of the Eurographics/ACM SIGGRAPH Symposium Geometry Processing, 2003, pp. 127–137. [18] X. Guo, X. Li, Y. Bao, X. Gu, H. Qin, Meshless thin-shell simulation based on global conformal parameterization, IEEE Trans. Vis. Comput. Graph. 12 (3) (2006) 375–385. [19] Y. He, M. Jin, X. Gu, H. Qin, A C 1 globally interpolatory spline of arbitrary topology, in: Proceedings of the 3rd IEEE Workshop on Variational, Geometric and Level Set Methods in Computer Vision, Lecture Notes in Computer Science, vol. 3752, 2005, pp. 295–306. [20] Y. He, K. Wang, H. Wang, X. Gu, H. Qin, Manifold T-spline, in: Proceedings of Geometric Modeling and Processing, 2006, pp. 409–422. [21] V. Jain, H. Zhang. Shape-based retrieval of articulated 3d models using spectral embedding, in: GMP, 2006, pp. 299–312. [22] D.L. James, C.D. Twigg, Skinning mesh animations, ACM Trans. Graph. 24 (3) (2005) 399–407. [23] M. Jin, Y. Wang, S.-T. Yau, X. Gu, Optimal global conformal surface parameterization, in: Proceedings of the Conference on Visualization’04, 2004, pp. 267–274. [24] S. Katz, A. Tal, Hierarchical mesh decomposition using fuzzy clustering and cuts, in: ACM SIGGRAPH 2003 Papers, 2003, pp. 954–961. [25] L. Kavan, J. Zˇára, Spherical blend skinning: a real-time deformation of articulated models, in: Proceedings of the 2005 Symposium on Interactive 3D Graphics and Games, 2005, pp. 9–16. [26] P.G. Kry, D.L. James, D.K. Pai, Eigenskin: real time large deformation character skinning in hardware, in: Proceedings of the 2002 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2002, pp. 153–159. [27] T. Kurihara, N. Miyata, Modeling deformable human hands from medical images, in: Proceedings of the 2004 ACM SIGGRAPH/ Eurographics Symposium on Computer animation, 2004, pp. 355– 363. [28] T.-Y. Lee, C.-H. Lin, H.-K. Chu, Y.-S. Wang, S.-W. Yen, C.-R. Tsai, Mesh pose-editing using examples, Comput. Animat. Virtual Worlds 18 (4–5) (2007) 235–245. [29] J.P. Lewis, M. Cordner, N. Fong, Pose space deformation: a unified approach to shape interpolation and skeleton-driven deformation, in: SIGGRAPH, 2000, pp. 165–172. [30] J.-M. Lien, J. Keyser, N.M. Amato, Simultaneous shape decomposition and skeletonization, in: Proceedings of the 2006 ACM Symposium on Solid and Physical Modeling, 2006, pp. 219–228. [31] W.-C. Ma, F.-C. Wu, M. Ouhyoung, Skeleton extraction of 3d objects with radial basis functions, in: Proceedings of the Shape Modeling International 2003, 2003, p. 207. [32] N. Magnenat-Thalmann, R. Laperrière, D. Thalmann, Joint-dependent local deformations for hand animation and object grasping, Proc. Graph. Interface 88 (1988) 26–33. [33] P. Merrell, Example-based model synthesis, in: Proceedings of the 2007 Symposium on I3D, 2007, pp. 105–112. [34] B. Merry, P. Marais, J. Gain, Animation space: a truly linear framework for character animation, ACM Trans. Graph. 25 (4) (2006) 1400–1423. [35] A. Mohr, M. Gleicher, Building efficient, accurate character skins from examples, in: ACM SIGGRAPH, 2003, pp. 562–568. [36] M. Mortara, G. Patanè, Affine-invariant skeleton of 3d shapes, in: Shape Modeling International, 2002. pp. 245–252. [37] J. Mukherjee, B. Chatterji, Thinning of 3d images using the safe point thinning algorithm (spta), Pattern Recogn. Lett. 10 (1989) 167–173. [38] X. Ni, M. Garland, J.C. Hart, Fair morse functions for extracting the topological structure of a surface mesh, ACM Trans. Graph. 23 (3) (2004) 613–622. [39] T. Oda, Y. Itoh, W. Nakai, K. Nomura, Y. Kitamura, F. Kishino, Interactive skeleton extraction for 3d animation using geodesic distances, in: ACM SIG, 2006, p. 9. [40] V. Pascucci, G. Scorzelli, P.-T. Bremer, A. Mascarenhas, Robust on-line computation of reeb graphs: simplicity and speed, in: ACM, 2007, p. 58.

62

Y. He et al. / Graphical Models 71 (2009) 49–62

[41] T. Pavlidis, Thinning algorithm for discrete binary images, j-CGIP 13 (2) (1980) 142–157. [42] G. Reeb, Sur les points singuliers d’une forme de pfaff completement intergrable ou d’une fonction numerique [on the singular points of a complete integral pfaff form or of a numerical function], Comptes Rendus Acad. Science Paris 222 (1946) 847–849. [43] D. Reniers, A. Telea, Skeleton-based hierarchical shape segmentation, in: Proceedings of the IEEE International Conference on Shape Modeling and Applications, Washington, DC, USA, 2007, IEEE Computer Society, 2007, pp. 179–188. [44] R.M. Rustamov, Laplace–Beltrami eigenfunctions for deformation invariant shape representation, in: Symposium on Geometry Processing, 2007, pp. 225–233. [45] S. Schaefer, C. Yuksel, Example-based skeleton extraction, in: Proceedings of the Fifth Eurographics Symposium on Geometry Processing, 2007, pp. 153–162. [46] T.W. Sederberg, S.R. Parry, Free-form deformation of solid geometric models, in: SIGGRAPH, 1986, pp. 151–160. [47] A. Sharf, M. Alexa, D. Cohen-Or, Context-based surface completion, ACM Trans. Graph. 23 (3) (2004) 878–887. [48] P.-P.J. Sloan, I. Charles, F. Rose, M.F. Cohen, Shape by example, in: Proceedings of the 2001 Symposium on Interactive 3D Graphics, 2001, pp. 135–143.

[49] R.W. Sumner, J. Popovic´, Deformation transfer for triangle meshes, in: ACM SIGGRAPH, 2004, pp. 399–405. [50] C. Theobalt, E. deAguiar, M.A. Magnor, H. Theisel, H.-P. Seidel, Marker-free kinematic skeleton estimation from sequences of volume data, in: Proceedings of the ACM Symposium on Virtual Reality Software and Technology, 2004, pp. 57–64. [51] J. Tierny, J.-P. Vandeborre, M. Daoudi, 3D Mesh skeleton extraction using topological and geometrical analyses, in: 14th Pacific Conference on Computer Graphics and Applications (Pacific Graphics 2006), 2006, pp. 85–94. [52] Y. Tong, P. Alliez, D. Cohen-Steiner, M. Desbrun, Designing quadrangulations with discrete harmonic forms, in: Symposium on Geometry Processing, 2006, pp. 201–210. [53] O. Weber, O. Sorkine, Y. Lipman, C. Gotsman, Context-aware skeletal shape deformation, Compu. Graph. Forum (Proceedings of Eurographics) 26 (3) (2007). [54] Z. Wood, H. Hoppe, M. Desbrun, P. Schröder, Removing excess topology from isosurfaces, ACM Trans. Graph. 23 (2) (2004) 190– 208. [55] X. Xiao, J.P. Lewis, S. Soon, N. Fong, T. Feng, A powell optimization approach for example-based skinning in a production animation environment, in: Computer Animation and Social Agents, 2006, pp. 141–150.