1 Introduction
Whether through laser range scanners, stereo reconstruction or structured light, methods and systems for 3D sensing and acquisition are now commonplace. However, these systems typically incur some amount of noise. Further, if we are interested in a particular object within data, it may be occluded by other objects. To extract the shape of an object of a particular class from such data it is often advantageous to use a statistical shape model to ensure that the extracted shape is valid for that object class. Another way of saying this, is that the shape is restricted to lie in a statistical shape space. In this paper, we review and describe, within a common mathematical framework, a wide variety of statistical shape spaces used for robust fitting to noisy and ambiguous data. We further perform a thorough theoretical analysis and experimental comparison of variants of such a model in which statistics are learned over global and local extents.
Such a statistical model must be learned from a database of consistently parametrized (i.e. registered) instances from the object class in question. One class of shape for which there exist a number of available databases is the shape of the human face, and this is the shape which we use for evaluation in this paper. Furthermore, the 3D shape of human faces is important to a wide variety of applications, ranging from telepresence to virtual avatar control to face recognition. However, we emphasize that the principles, models and algorithms discussed in this paper are applicable to any class of shapes for which a database containing parametrized data is available.
The main reason to use a statistical shape model to fit to data, instead of a template fitting with a nonrigid iterative closest point (ICP) approach and some kind of regularization constraint, is that by learning a statistical model for a class of shapes, we can significantly reduce the search space, which results in the ability to reconstruct the underlying shape in the presence of severe noise or occlusions. These and other types of ambiguities are often present in realworld data captured in uncontrolled environments.
The purpose of this paper is to provide researchers, engineers and endapplication developers interested in employing existing statistical models or developing new ones with a targeted, focused review that includes a thorough analysis and comparison of the costs and benefits of two statistical models for the task of model fitting to noisy, corrupted or incomplete 3D data. The task of fitting a 3D shape model to ambiguous data is important for many applications, such as recognition tasks (identity and expression recognition in the case of faces, sex recognition and gait analysis in the case of human bodies), telepresence, virtual avatar control, segmenting and extracting organ shapes from medical images for surgical planning or diagnosis, and so on. We specifically maintain a targeted scope to provide an indepth analysis of the behaviors of various statistical models. Toward this end, we analyze different statistical models from those that model global geometry variation to those that model local geometry variation. We experimentally compare a purely global model to a purely local model, thereby providing a quantitative comparison of two models on either end of the spectrum of commonly used statistical models. Thus, the contribution of this paper is threefold:

An analysis of the theoretical properties, within a common mathematical framework, of a wide variety of seemingly very different statistical shape spaces.

A quantitative and qualitative analysis, and extensive comparative evaluation of the practical performance of both global and local statistical shape models.

We publish the learned statistical models and code to use them [1], thus allowing others to try them out and potentially fit them to other input modalities and for any application.
We begin by reviewing the statistical shape spaces, or statistical shape models, that have been used for human face shapes, human body shapes, and medical data (Section 2). We give a common mathematical framework in which these models relate to one another and discuss the key differences between them. We then give a mathematical description of the process of learning or training a statistical shape space, and an indepth focus on two specific models for human faces that represent extremes in the range of models (Section 3). This is followed by a mathematical description of the process of fitting a statistical model to ambiguous, partial or corrupted data via energy minimization (Section 4). Again, this is accompanied by an indepth analysis of the process for the two shape spaces. We then provide an extensive experimental comparison of the two models (Section 5), followed by a discussion of the implications of the results (Sections 6 and 7).
2 Statistical Shape Spaces
To extract an object from noisy input data, it helps to have a small basis in which the shape of the object can be represented. Traditionally, such bases have been generated for specific classes of models by artists. An example of such an artistgenerated basis are blendshape models (as for instance used by Li et al. [2]
), which can be used to encode a face performing different expressions. This way of generating a basis requires expertise about the possible deformations, and modeling the shapes is tedious. More recently, machine learning has been used to find a small basis from a set of training shapes using statistical analysis. This gives an easy and fully automatic way to find a small set of basis functions.
In this vein, we formulate statistical shape analysis of a given object class as the task of finding a statistical shape space that efficiently and informatively represents the shape of objects of that class. We define a statistical shape space as a shape space
equipped with a probability distribution, or
prior, measuring how likely it is that an object of the given class would have a particular parametric representation in the shape space. The shape space itself is defined by the set of coefficients obtained by projecting the shapes onto the set of basis functions. (We use the terms basis functions and basis vectors somewhat interchangeably in the following; strictly speaking a basis function is only relevant for continous surfaces, and in practice basis vectors are used for discrete data.)
Thus, we focus on statistical shape analysis as a generative technique. A surface containing vertices in is represented by shape parameters or coefficients, which form a vector . A generator function
(1) 
generates from these shape parameters a surface representation (either mesh or point cloud) of vertices. These shape parameters, and by extension the surface, can be fit to input data of varying modalities (3D point clouds, 3D voxel images, 2.5D depth, 2D images, sparse measurements, etc.), so long as there is a way to measure the distance between the surface and the data, or the quality of the fitting.
As we see in the following review, by far the most common form of statistical analysis used for shapes is principle component analysis (PCA), which seeks a basis in which variance of the training data is maximized. The resulting basis vectors are the directions of greatest variation within the training data. Projecting the training samples onto this basis results in a diagonal sample covariance matrix. If the underlying distribution of the data is assumed to be multidimensional Gaussian, then this corresponds to the maximum likelihood estimate of the parameters of the density function. As a result, the resulting shape space is often equipped with a Gaussian prior. If this assumption does not hold, then a Gaussian prior may be arbitrarily far from the true prior, and choosing the correct prior may be challenging. A mathematical description of how to learn a statistical shape space using PCA is given in Section
3.2.The core aspect that varies between statistical shape spaces is the space in which PCA is performed. Careful selection of the space in which to perform PCA can lead to significantly better statistical properties. For many models, this amounts to a change of basis, followed by separate analyses in different subspaces of the transformed space. The resulting basis of the learned space is then formed by composing the decomposition basis with the subspace PCA basis (directions of greatest variation within the subspace). For others, this amounts to a nonlinear transformation followed by a global analysis.
This section reviews work on performing statistical shape analysis of 3D data for image processing applications. The categorization of the statistical models in different application domains is summarized in Table 1, where models and methods are grouped by the type of data they were applied to and by the extent of the basis functions used. The table further contains information on whether the models were designed to analyze articulation variations separately from shape variations (here, articulation can refer to facial expression, body posture, or the pose of bones before and after an operation).
In this study, we focus on shape variations over a sample from a population, and hence in Sections 3.2, 3.3 and we analyze and compare statistical models for faces without expression variations. In Section 5 we perform an extensive experimental comparison of the models. This allows us to better examine the differences between the statistical models themselves with respect to the modelfitting task.
The first step to performing shape analysis is to acquire and register a set of training shapes that capture the shape variability that is of interest for a particular application. Subsequently, statistical shape analysis is performed on the registered training shapes: the shapes are projected onto a basis of choice and a probability distribution is fitted to the resulting coefficients to obtain a prior distribution for the shapes of interest. Without correspondence information, this statistical analysis is not possible. However, as indicated in the last column of Table 1, a few methods simultaneously compute a parametrization of a population of shapes while building a statistical model.
Computing correspondences between a population of shapes is a challenging problem, and a detailed discussion about possible approaches is beyond the scope of this work. We refer the reader to recent surveys [3, 4] for more information. However, we emphasize that the quality of the registration greatly affects the quality of the resulting statistical models, and by using a highquality registration in this study, computed as discussed in Section 5.1, we are able to better analyze the properties of models themselves rather than the effects of gross misregistration.
In computer vision, statistical 3D shape models are commonly used to infer the threedimensional shape of an object from images, mostly for the purpose of image manipulation. While recently, different classes of shapes have been considered
[5, 6], shape models of human faces and human bodies are of special interest due to their immense applicability in human–machine interaction. In medical image analysis, statistical shape models are commonly used to segment medical images and to find correspondences and abnormalities of anatomical shapes. In the following, we review statistical shape spaces used to analyze human faces (Section 2.1), human bodies (Section 2.2) and medical data (Section 2.3).Type of data  Influence of  Methods  Articulation variation  Simultaneous 
basis functions  analyzed separately  parametrization  
Faces  
global  Morphable model (PCA) [7, 8, 9, 10]  
global  Morphable model (PCA) [11]  
global  Multilinear model [12, 13, 14, 15]  
partbased  Partbased model [16, 17, 18]  
partbased  Partbased model [19]  
localized detail  Hierarchical pyramids [20]  
local  Local wavelet model [21]  
Bodies  
global  PCA model [22, 23, 24, 25, 26]  
global  SCAPE model [27, 28, 29, 30, 31, 32]  
global  SCAPE model [33]  
global  Rotationinvariant encoding [34, 35]  
global  Multilinear model [36]  
global  Postureinvariant model [37]  
partbased  Segmented PCA model [38]  
partbased  Partbased multilinear model [39]  
Medical Data  
global  Active shape model (PCA) [40, 41]  
global  Active shape model (PCA) [42]  
global  PGA model [43]  
partbased  Partbased model [44]  
partbased  Partbased multilinear model [45]  
local  Local wavelet model [46, 47, 48, 49, 50, 51] 
2.1 Human Face Shapes
We start our review by summarizing the use of statistical shape models of human faces. Blanz and Vetter [7] proposed the first statistical shape model for 3D models of human faces. The model, called morphable model, captures both 3D shape and texture information and can be used to predict a 3D face shape from a single input image. Texture is an important cue for human faces, and can greatly help with fitting directly to images. However, in this paper, we focus on statistical analysis of shape for clarity. A parametrized database of 3D face scans, mostly in neutral expression, is used to learn the statistical model using standard PCA. Given an input image in neutral expression, the learned shape space is searched to find the textured shape that best explains the input image using an optimization technique. This successful approach resulted in multiple followup works [8, 9]. Amberg et al. [11] extended the morphable model to model expression variation by computing offsets from the neutral pose and performing PCA over these offset surfaces. With this method, they are able to include expression variation while maintaining a single linear model.
Blanz and Vetter [7] experimented with manually segmenting the morphable model into four regions. A morphable model is then learned for each region and the regions are fitted to the data independently and merged in a postprocessing step. This partbased model was shown to lead to a higher data accuracy than the global morphable model [7]. As this approach is suitable to obtain good fitting results in localized regions, it has been used in multiple followup works [16, 19, 17]. All of these methods proceed by manually segmenting the faces, by learning an independent morphable model for each part, and by fitting the parts to data independently. In this case, a basis arises out of the set of indicator functions on the segments, which are then composed with the union of PCA bases on the segments. As a result, the areas at the boundaries of the individual parts are not necessarily continuous, and a postprocessing step is used to merge the fitted patches. Smet and Van Gool [18] proposed a similar method where the segmentation is found automatically by clustering the vertices based on features derived from their displacements over the training set. To address the potential discontinuities at the boundaries of the segments, they smoothly weight the segments to obtain regionalized basis functions for the training data. To estimate the optimal weights, a complex iterative algorithm is used. In each iteration, the training set is randomly partitioned into two disjoint subsets, where one set is used as test set to be reconstructed from the model learned on the other. This is followed by a twostep procedure. First, optimal reconstruction coefficients are estimated by weighted leastsquares. Second, the optimal weights are estimated by solving many independent linear systems. Since for highresolution training data this iterative algorithm converges slowly, they employ a coarsetofine approach to accelerate convergence.
Brunton et al. [21] used a statistical analysis based on wavelet models to learn many localized independent prior distributions at multiple scales for the 3D shape of human faces. This allows to capture and combine localized shape variations in different areas of the face in a multiresolution framework independently. They used this information to predict the 3D face shape from point clouds generated by stereo reconstruction. In this model, a basis is formed by composing the wavelet basis with the local PCA bases. As in the work by Blanz and Vetter, all faces are assumed to have a neutral facial expression. While PCA models such as the morphable model and partbased PCA models are commonly used in computer vision to model human faces, localized prior distributions based on wavelet models are only beginning to influence this field. In the following, we will refer to this method as local wavelet model.
Golovinskiy et al. [20] proposed a statistical model based on hierarchical pyramids to synthesize geometric facial details. This statistical model, which allows for spatially varying geometric detail across the face, models the difference between a smooth face and a highresolution face with geometric details, such as wrinkles.
To allow for varying facial expressions, Vlasic et al. [12]
used a statistical method based on tensor algebra called the
multilinear model to analyze a set of faces captured of subjects performing a variety of facial expressions. This approach can be viewed as an extension of morphable models to higher dimensions. It is formed by taking the Cartesian product of two linear statistical shape spaces that capture variations according to two different attributes (i.e. identity and expression). The resulting basis is then the Cartesian product of the two bases associated with each attribute. This model allows for the prediction of a 3D face shape in multiple possible expressions from a single photograph.Yang et al. [10] exchanged the expression of a face in a single image based on a different input image of the same subject. For this application, they built multiple PCA spaces (one per expression) and combined these spaces for their application. A followup paper [14] used a multilinear model to enhance or dampen expressions in videos.
More recently, much work has focused on extracting a set of frames of threedimensional face shapes from a video stream showing a face. The output of this type of algorithm is a fourdimensional sequence showing the threedimensional face shape in motion. Dale et al. [13] extended the method by Vlasic et al. to compute such a fourdimensional sequence. This information is then used to exchange faces in video sequences; either keeping the sequences of expressions and replacing the identity of the subject, or keeping the identity and replacing the sequence of expressions. Bolkart and Wuhrer [15] learn a multilinear model and fit it to sequences of 3D faces performing various expressions, producing a 4D parametrization, which can be used to animate a static scan with a specified expression.
Another avenue of recent work is to track twodimensional range images over time. These images can be captured using depth sensors, such as the Kinect sensor. Weise et al. [52] proposed a method to track a threedimensional face model over time using prior information on the deformation model, and to use the tracked model to drive the animation of a virtual character in realtime. Unlike the other methods discussed in this section, this method learns a statistical prior that is subjectspecific. In this paper, we do not consider subjectspecific priors.
2.2 Human Body Shapes
Allen et al. [22] proposed a statistical model for human body shapes acquired in a similar posture that is similar to the morphable face model introduced by Blanz and Vetter. One main difference is that Allen et al. only learn information about the 3D shape and not about texture. We denote this model by global PCA model in the following. This model has been used to predict a 3D human body shape in a similar posture from one or more images [23, 24, 25] and from measurements [26]. Xi et al. [38] used a partbased PCA model to infer body shapes from silhouettes. Wuhrer et al. [37]
performed PCA on a local shape representation based on the Laplace operator to analyze human body shape variations independently of posture changes. This is an example of a nonlinear transform followed by a global PCA in the transformed space, where the transformed space has been chosen to be invariant to local rigid deformations, making the resulting statistical analysis posture invariant.
To allow for posture variation, Anguelov et al. [27] proposed the SCAPE model. This model learns a PCA shape space for body shape variations using a database containing multiple subjects in a similar posture. Furthermore, the model learns a mapping from posture parameters (based on a skeleton) to shape changes using a database containing one subject in multiple poses. The model then combines the two variations using the assumption that body shape and posture are decorrelated. Since the SCAPE model successfully models human body shape and posture, it has been used to predict a 3D body shape in arbitrary posture from a single image [28]. Furthermore, this model can be used to predict a 3D human body shape in arbitrary posture based on a set of input images of a dressed person [29]. Such a 3D prediction can then be used to modify the input image [30]. Just like in the case of human faces, more recently, much work has focused on finding a fourdimensional sequence of threedimensional human body shapes in motion from a video sequence. Jain et al. [31] used this to modify human body shapes in video sequences. Weiss et al. [32] proposed to use the SCAPE model to compute a 3D body scan from noisy Kinect data.
A recent approach by Hirshberg et al. [33] treats the problems of learning a SCAPE model and computing pointtopoint correspondences between a set of training data simultaneously by solving a single variational problem.
A different avenue to allow for posture variation is to model shape and posture changes as correlated. This assumption is relevant, since the difference of the human body shape of the same subject in different postures depends on the body shape, e.g. on how muscular the subject is. Hasler et al. [34] proposed a shape space that jointly captures shape and posture variations by performing PCA on a rotationinvariant encoding of the shapes. This shape space is then used to predict the body shape of a dressed subject [35].
An alternative for correlating human shape and posture variations is to use a multilinear model (as Vlasic et al. [12] do for face shapes). This avenue was explored by Hasler et al. [36] for the application of predicting 3D body shape and posture from an image. Chen et al. [39] proposed a partbased multilinear model using a manual segmentation into body parts.
2.3 Medical Data
In medical imaging, one is especially interested in a body part, such as an organ or part thereof. Tasks of interest include finding the shape of interest in a medical image. This decomposition of the image is often called segmentation. To solve this task, Cootes et al. [40] propose the use of a statistical prior called active shape model, which is similar to the morphable model for faces. The active shape model learns the distribution of a set of registered and aligned training shapes using PCA, and uses this prior information to segment a given medical image. This model is commonly used for image segmentation, see Cootes and Taylor [41] and references therein.
The quality of the registration of the training shapes directly influences the quality of the statistical shape model. Davies et al. [42] used this observation to derive an approach that jointly optimizes the registration and a shape model built using PCA. The main idea of this approach is that a good shape model should have a small informationtheoretic description length. Hence, to simultaneously parameterize the shapes and build a statistical model, the approach by Davies et al. solves a variational problem that aims to minimize the description length of the PCA model.
PCA assumes that the data can be approximated well using a linear model. However, in medical imaging, many data sets form a nonlinear manifold in a highdimensional space. To capture the structure of this highdimensional manifold, Fletcher et al. [43] generalized PCA to this setting. This approach called principal geodesic analysis (PGA) considers geodesic distances between shapes measured along the highdimensional manifold instead of Euclidean distances.
One problem with active shape models is that they capture global shape variations. In medical imaging, one is often interested in detecting localized shape anomalies, as these can give insights in whether or not a specific organ is affected by a disease, for instance. In an active shape model, such local variations may be distributed over several principal components, and they may be controlled by principal components that capture a small percentage of the overall shape variability.
To remedy this, partbased models have been proposed for medical imaging. Toews et al. [44] proposed a partbased model that captures the statistical variations of geometry and color. This model is used to analyze MRI volumes of the brain. Lecron et al. [45] used a different partbased model to analyze variations of the spine both within and across subjects. Each vertebra is considered one part and analyzed using a multilinear statistical model.
Partbased models work well for applications where a segmentation into parts is meaningful. However, in many medical imaging tasks, a natural segmentation is difficult to define. To address this problem, Davatzikos et al. [46] used statistical analysis based on local wavelet models to learn a localized prior distribution of contours in images. Here, the localized regions with large variations are found automatically and do not need to be predefined. Nain et al. [47] extended this technique to use wavelets to perform a statistical analysis of threedimensional shapes. Shape priors based on different types of wavelet models have been used to segment medical images [48, 49, 50]. Yu et al. [51] show that statistical wavelet models can be used to analyze cortical folding patterns, which is a challenging task.
2.4 Comparison
In the following, we provide an analytical comparative evaluation of a global PCA model and a local wavelet model. By doing so, we compare two methods on either end of the spectrum of the commonly used statistical shape models. It is already known [7] that manually segmented partbased decompositions improve fitting accuracy for human face shapes. A thorough numerical comparison has not been performed for automatically and objectclass nonspecific decompositions such as wavelets. We expect that, like partbased models, the local wavelet model will obtain locally more accurate results than the global model, since the basis functions have limited extent. However, since the partdecomposition for statistical models is to our knowledge always shapespecific, the local model can more readily be adapted to various application areas. In particular, when segmenting volumetric medical data for a particular organ or part of the brain, it may be impossible to identify distinct parts on which partspecific statistical models can be learned. However, one may still want localized fitting for these shapes.
Applications of statistical model fitting are wideranging and require very different evaluation criteria. To evaluate them all would be impractical, and to evaluate one or two would unduly reduce the scope of this study. Our comparison (Section 5) therefore uses as a test metaapplication the task of model fitting to point cloud data of human faces, and from the results we give an intuition of which model would be most appropriate for which end application (Section 6). All of the figures and tables in the following sections depend on the specific training data used, but they nonetheless illustrate how the geometric information captured in the statistical models can be evaluated and visualized. This comparison can provide a guide of how the locality of the basis functions of a model influence its performance for model fitting. Reconstruction of 3D shape from silhouettes or images, or segmentation of volumetric data, can also be viewed as a model fitting, with additional or different ambiguities, and therefore we infer that we can reasonably expect the statistical models to behave similarly in these scenarios.
Furthermore, we allow those developing end applications to test the performance of both models in their specific scenario by publishing the statistical models and code to use them online [1]. The data and code provided will also allow others to derive fitting energies and code for different input modalities.
3 Learning a Statistical Shape Space
For learning a statistical shape space, or training a statistical shape model, we assume we are given training shapes in full correspondence. Figure 1
gives an overview of the model training. The key difference between the two models discussed in detail here is the basis in which a prior probability distribution is fit to the training data.
We prealign the data to remove rotation, translation, and uniform scale differences using generalized Procrustes analysis (GPA) [53]. Note that by removing uniform scale differences, we only consider shape differences and not size differences of the models. For datafitting this is desirable due to different measurement units used by different acquisition systems, but in general this is application dependent [53]. GPA iteratively aligns each model to the mean shape and recomputes the mean. Removing transformations that are not of interest using GPA is an important preprocessing step that yields better statistical models.
Learning a statistical shape space requires determining the basis functions or vectors, which define the shape space, and fitting a probability distribution to the resulting shape space coefficients from the training set. PCAbased methods perform these two steps simultaneously, selecting as a basis the directions of greatest variations in the data and computing a diagonal sample covariance matrix for the data projected onto this basis, which corresponds to a maximum likelihood estimate of a multidimensional Gaussian distribution. Partbased and waveletdomain methods decompose the shapes into a localized basis before proceeding with PCA. The result is a basis consisting of the localized basis functions composed with the localized principal components. The learned prior is the product of the localized multidimensional Gaussian distributions. Note that the assumption of a Gaussian density function is only introduced when equipping the shape space with a Gaussian prior.
With this framework, and by restricting ourselves to linear shape spaces for the purposes of this study, the generator function Eq. (1), can be written as a combination of the basis functions
(2) 
where is the mean shape computed over the training set, is a matrix, are its columns, and as before is a vector of shape parameters. It is precisely the choice of that determines the properties of the shape space, and determines the prior distribution learned from the training samples.
3.1 Evaluating a Statistical Model
If a statistical model with a small number of basis functions is fitted to a face, the result contains little shape detail, because the model only represents a small proportion of the variability of the training data. Keeping a large number of basis functions may cause the model to overfit to the training data. That is, the learned space may contain a bias towards shapes present in the training set. Overfitting can occur when the model is underconstrained by the training samples. This may occur if the model itself has too many degrees of freedom, or when it is learned in a very highdimensional space, such as when the model is computed directly from very highdimensional training samples.
To pick a number of basis functions that preserves a high amount of variability yet does not overfit the training data, we use the following three error measures similar to compactness, generalization, and specificity [54]. We use a slight modification of the original error measures to obtain results that are independent of the size of the training data.
Compactness measures how much variability of the training data is explained by the learned statistical model. That is, we want to measure what fraction of the total variability of the training data is captured by model parameters. This provides a measure of how well a given number of parameters explains the training data.
Generalization
measures the ability of the model to represent data, which are not part of the training set. To calculate this measure, we learn a PCA model on a subset of the training data, where one subject is excluded. The excluded subject is projected to the PCA space, reconstructed, and the distance between the source and the reconstruction is measured. To measure the distance between two faces, we use the average Euclidean vertex distance computed between all corresponding vertices. We perform this measurement for all subjects. The mean and standard deviation are then considered.
Specificity measures the similarity between reconstructions from the statistical model and the training data. This estimates the plausibility of a random face represented using the learned shape space. To calculate specificity we choose a set of random points sampled from the probability distribution of the learned statistical shape space. For each of these points we reconstruct the shape using Eq. (2) and compute the distance to the closest face in the training data. The distance between two faces is computed as above. The mean and standard deviation for the random sample are then considered.
3.2 Global PCA
Principal component analysis aims to reduce the complexity of a set of data. Due to its simplicity it is widely used for shape analysis. PCA is a linear transformation of a set of vectors from to with . A vector is expressed by the scalar weights in a dimensional subspace, spanned by the orthogonal vectors , by
(3) 
For each parameterized shape of the training set we have one vector that contains an ordered coordinate set of all points of the th training shape. The vectors
are the eigenvectors of the data covariance matrix
(4) 
where is the mean of the training data. The eigenvectors
are ordered with respect to the nonincreasing corresponding eigenvalues
. The eigenvalues measure the variability captured by the th principal component. More specifically, captures of the variability of the training data. The rank of the data covariance matrix is at most and therefore the number of distinct nonzero eigenvalues and hence, the number or principal components, is at most .Thus, we get our basis directly from the data via the principal components: in matrix form has columns , for , where . Note that every basis vector has global support in general: all elements are in general nonzero, and every vertex is influenced.
Properties The global shape space represents the highdimensional differences of the training faces in a lowdimensional shape space that is spanned by the corresponding eigenvectors of the largest eigenvalues of the data covariance matrix.
Figure 2 visualizes the variations along two principal components in the range of to , where denotes the standard deviation of the th principal component.
The amount of details that can be expressed by the global statistical model is limited by the details present in the training data. It would be useful in some applications to increase the variability of the training data by increasing the mesh resolution of the training data by inserting new vertices as linear combinations of existing vertices. Unfortunately, this is not possible using this global approach. If a vertex expressed as a fixed linear combination of existing vertices is inserted to each training mesh, the corresponding additional vertex in the fitted surface is identical to the corresponding fixed linear combination. Therefore, if an additional vertex is chosen to be placed on a triangle, the corresponding additional vertex is located on the corresponding triangle of the fitted result. Hence, using fixed linear combinations to add points to the surface of the model and fitting this extended surface to a target face leads to the same result as fitting the original model to the target face and adding the points into the resulting surface using the fixed linear combinations. This is a key difference to the local model reviewed in Section 3.3.
Computing the compactness statistical measure is straightforward for the global PCA model. For principal components, compactness is defined as
(5) 
where is the th eigenvalue of the data covariance matrix. Computation of generalization and specificity are also straightforward.
3.3 Wavelet PCA
The core principle behind wavelet transforms is to project sampled data onto a set of basis functions that are localized in space and frequency. In our context, we use the wavelet basis as a prefix, and extract data driven basis for individual wavelet coefficients using PCA. Secondgeneration or lifting wavelets [55] are computed in time linear in the number of samples in the original signal using local lifting operations and subsampling at each scale. The samples are partitioned into maximally correlated subsets, e.g.odd and even samples for signals on 1D domains. One subset is then used to predict the other, and the residual of this prediction is then called the detail, or wavelet, coefficients. The detail coefficients are then used to update the other subset, giving approximation, or scaling, coefficients. The process is repeated on the scaling coefficients. We refer the reader to Sweldens [55] for a more thorough explanation.
While wavelets were originally defined regularly sampled Euclidean domains [56], spherical wavelets [57] are defined on subdivision surfaces, often topological spheres. A commonly used wavelet basis is a biorthogonal generalized Bspline basis [58] that uses the CatmullClark subdivision scheme and has been applied in multiple application domains [49, 21]
. The prediction and update operators are Bspline interpolations from the neighboring vertices. This scheme is stable for linear and cubic Bsplines; our comparison uses the linear basis as we found that this produced satisfactory results.
Aside from wavelet coefficients, one might use any localized basis as a prefix to PCA or other datadriven bases. However, wavelet bases have the important property that they decorrelate the data, meaning the resulting coefficients can be analyzed separately. Since individual coefficients are of much lower dimension than the whole surface, this greatly reduces the risk of overfitting, given the same number of training samples .
Performing PCA over the whole set of wavelet coefficients would result in the same principal components as the global model, because the wavelet transform is a linear transform, and PCA essentially just rotates the data so that the coordinate axes align with the directions of greatest variations. Instead, this method performs PCA locally on each coefficient, which is a 3D vector quantity, over the database.
First, let us denote the mean of each wavelet coefficient over the database by
(6) 
where indexes the coefficients.
While we can perform statistical analysis on each independently of other values of , we must consider their three components together. Each is a 3D vector representing either the scale (absolute value) or the detail (relative value) of the shape at a particular frequency and spatial location. However, the coordinate axes in general do not correspond to the directions of greatest variation in the database. Therefore, we perform PCA on each set of coefficient vectors, to obtain 3D vectors that represent the position along the directions of greatest variation, and matrices that transform these coordinates to our original world coordinate system, as in
(7) 
where we write and to denote the components of these vectors. Applying the transform to the data diagonalizes the covariance matrix, thus making each component independent.
The reconstruction of the face shape from the model is then given by the inverse wavelet transform
(8) 
where is the vertex index in the reconstructed surface, is the level of wavelet coefficient, is the number of levels used, is the set of scaling functions at level zero, is the set of wavelet functions at level , and are the coefficient indices, is the scaling function at the coarsest level centered on location , and is the wavelet function at level and location . While the transform is expressed here in terms of basis functions, it is computed using lifting operators, which amount to weighted averages of a vertex’s local neighborhood, and it can be expressed as a matrix multiplication. The basis functions themselves, and , are Bspline approximations to Gaussian and Mexicanhat functions. For more details, see Bertram et al. [58].
We can now construct the combined basis as follows. Observing that the inverse wavelet transform, Eq. (8), is a linear operator on the vector of concatenated wavelet coefficients s, we can write it as a matrix . Thus, we have
(9) 
and from Eq. (7) we have
(10) 
where is the concatenation of the coefficient means , and is a blockdiagonal matrix with the matrices on the diagonal. Therefore, because , we have
(11) 
as our combined basis, and the dimensionality of our shape space is . Note that has full rank.
To use a spherical wavelet basis to represent shape, it must be a subdivision surface, in our case CatmullClark subdivision hierarchy. For training, the surfaces in the database are typically stored as triangle meshes without subdivision structure. Thus, we must resample the surfaces with the proper structure. The subdivision scheme uses quadrilateral elements, although it can handle extraordinary vertices. We resample the triangle meshes using the custom, yet straightforward, technique by Brunton et al. [21] tailored to the fact that we are dealing with faces, which are topologically like a disc. In principle, however, any quadremeshing technique can be used.
Properties The local model has the benefit that it avoids overfitting, and as a consequence we can keep all variability present in the training set. Intuitively, the local surface properties of any given surface point are not likely to be specific to one set of faces or another. Whereas for the global model a bias in the training set, overrepresentation of one sex or a particular ethnicity or age range, can cause the lesser principal components to be highly specialized to that set, the geometry of a local surface patch is likely to be less dependent on the training data. The consequence is a somewhat unexpected behavior: by training and combining many lowdimensional models, which due to the limited flexibility of the training space () have reduced sensitivity to bias in the training set, we get a final model with much greater flexibility, because truncation becomes unnecessary.
Figure 3 visualizes the mean shape colorcoded with the magnitude of the shape variability for four levels of the wavelet subdivision, which corresponds to the localized shape variations at different scales. At finer scales, the variation quickly localizes around major facial features and reduces in magnitude.
level 0  level 1  level 2  level 3 
The dimensionality of the local model is statistically more favorable. If, as is usually the case, the number of vertices is much greater than the number of training examples, , then the global model has problems of fitting to the particularities of the training set. In the local model, many independent statistical priors are learned, each with dimension 3. We have many more training examples than that. The independence of the local priors further allows an exhaustive search of the parameter space. Thus, we have no danger of getting trapped in local minima.
The drawback of these properties, in particular of retaining all the variability of the training data, is that the local model is a much higherdimensional representation than the global model. Thus, the dimensionality of the local model can be computationally much less favorable. There is, however, a tradeoff that can be made by using the wavelet basis progressively, and working with the localized priors only up to a certain level.
As the lifting operations of the wavelet transform amount to local weighted averages of vertex coordinates, the transform can be expressed as a matrix multiplication, if the surface is expressed as a vector containing the vertex coordinates. Because the transform is biorthogonal, this matrix is square and has full rank. In contrast, the global PCA basis , has rank . As discussed in Section 3.2, resampling the surface at linear combinations of vertex coordinates (eg., within a triangle), does not increase the rank of the transform. However, because has full rank, we can obtain more detail by linearly upsampling the training surfaces. This means that we can resample the training shapes at high resolution to obtain a statistical model that captures fine shape detail.
The statistical measures generalization and specificity can be computed in the same manner as for the global PCA model. While there is no simple formula for the compactness of the wavelet model, if we retain all shape parameters, and hence all variation, the compactness measure is fixed at . In this case, the generalization measure is .
4 Fitting a Statistical Shape Model
In this section we give an overview of the generic model fitting approach shown in Figure 4. This assumes that a statistical shape model has been learned, for instance using one of the approaches given in Sections 3.2 and 3.3. As can be seen from Figure 4, the process begins with an initial alignment using either landmarks or feature points, described in Section 4.1, regardless of the statistical shape space. Subsequently, we minimize an energy function with respect to the model parameters (Section 4.2), which are specific to the statistical shape space, and may affect the minimization strategy (Sections 4.3 and 4.4).
(a)  (b) 
4.1 Initial Alignment
To fit a statistical shape model to an input data set, we first need to align the input data and the statistical shape model to be in the same global coordinate system. Since we consider only shape differences in the training data, the initial alignment aims to find the rotation, translation, and uniform scaling that best aligns the statistical shape model with the input data.
To compute such an initial alignment, corresponding landmarks are commonly used. These landmarks can be manually located on the mean shape of the aligned training database once. On the input data, the landmarks can be predicted in a fully automatic way [59, 60]. However, since we use a test database that contains a set of landmarks, we choose to use a subset of these landmarks (the ones shown in red in Figure 5a) to compute an initial alignment. This approach removes a potential source of fitting error due to landmark prediction inaccuracies.
Another commonly used way to rigidly align two shapes is to use automatically detected features. We test a method of this flavor in our experiments. The method we use proceeds by finding corresponding features on the mean face and the input scan using spin images [61] and by performing random sample consensus [62]. This fully automatic method is expected to lead to less accurate alignments than the use of the given landmarks.
4.2 Energy Minimization in Shape Space
Our goal is to fit the statistical shape model to the input data as closely as possible while staying in the learned shape space. To fit our model to data, we minimize an energy function that amounts to the sum of squared distances between each model vertex and its nearest neighbor in the input point cloud. For our experiments, we use the following commonly used basic energy to pull the model towards the data
(12) 
where is vertex of (see Eq. (1) and (2)), is a point in the input point cloud, returns the index of the nearest neighbor in of , and
is a truncation threshold to add robustness against outliers. We compute nearest neighbors with a kd tree using the implementation in ANN
[63].When fitting a statistical shape model to data, the space of possible solutions should only contain likely shapes, thus ensuring that only plausible results are possible. A common and intuitive approach is to use the (negative logarithm of the) learned prior distribution as an energy term. In the case of PCA, it is common to assume a multidimensional Gaussian centered on the mean shape. In terms of energy minimization, this amounts to placing a soft constraint that the solution should be close to the mean. By design, however, this introduces a bias into the optimization, and results using this technique tend to lose distinctiveness and look similar to the mean.
Patel and Smith [64] proposed an alternative prior that is aimed at maintaining the distinctiveness of the models. They model the shape space as a manifold that is at a constant Mahalanobis distance from the mean. This is based on the observation that the squared Mahalanobis distances from the mean of a set of
dimensional normally distributed vectors form a
distribution with expected value equal to . Hence, Patel and Smith restrict the shapes to be on the hyperellipsoid at Mahalanobis distance from the mean in order to preserve shape distinctiveness. While this approach models distinctiveness using the expected Mahalanobis distance from the mean, it does not consider the normal distributions along each dimension of the shape space. That is, the modeled shape space may contain highly unlikely shapes along the directions of the principal components, as can be seen for the shape at the intersection of the hyperellipsoid shown in red and the axis in Figure 6.To simultaneously avoid meanshape bias and highly unlikely shapes, in our experiments we constrain the shape to lie within a region of the shape space where the prior probability is sufficiently high, i.e. the shape is sufficiently likely. Ideally, in the case of PCA models, one would impose ellipsoidal constraints corresponding to a probability isosurface of a given value. However, from an optimization standpoint, it is much easier and more efficient to impose linear constraints. Hence, we constrain the shape to lie within the hyperbox of about the mean shape, where is the standard deviation of the training data along dimension of the learned statistical space, and is a parameter controlling the amount of deviation allowed. Figure 6 shows a twodimensional plot of this hyperbox. We demonstrate that in this shape space, by minimizing an energy function with only the data term given in Eq. (12), we can maintain distinctiveness, while avoiding unrealistic shapes. This is equivalent to a prior probability of the form
(13) 
where
(14) 
if we assume the shape parameters are centered (mean subtracted). We call this a hyperbox prior.
A common postprocess to fitting the parameters of the statistical models, is to then leave the statistical shape space and perform a finefitting of the vertex positions directly, similar to a template fitting method. We deliberately do not do this for two reasons. First, such a step is most often necessary when the learned shape space is not sufficiently generalizable to express novel shapes. This can occur due to insufficient or poor training data. However, as detailed in Section 5, we train from clean data containing a good sampling of both sexes and different ethnicities with a highquality registration. Thus, our learned models are of highquality.
Second, we wish to study the properties of the statistical shape spaces themselves, and leaving the space in a postprocess would inevitably introduce additional uncertainty in analyzing the fitting results in Section 5. This is particularly true when we evaluate the fitting in the presence of occlusions.
4.3 Global PCA
Energy Minimization While the energy in Eq. (12) is not strictly differentiable at all points, it is continuous, and the number of points where it is not differentiable is small. Hence, we can minimize it using a bounded QuasiNewton method [65]. The coordinate bounds on the parameters enforce the condition given in Eq. (14). This minimizer gives superlinear convergence rates without the need for an explicit inverse Hessian. Note that this optimization technique does not guarantee to find the global optimum of the energy function.
Computational Complexity Let denote the number of iterations required for the minimization to converge, and let denote the number of data points in the target shape. The complexity of building a kd tree of points in is , and a single nearest neighbor search takes time [66]. Given the nearest neighbor indices for all points, a single evaluation of the energy given in Eq. (12) takes time, and a single evaluation of its gradient takes time. Thus, the overall time complexity of fitting the global PCA model to a data set is . For this model, each training shape contains vertices.
4.4 Wavelet PCA
Energy Minimization We minimize the energy in Eq. (12) using a global search of each parameter. That is, we sample uniformly within the range given by the hyperbox prior given in Eq. (14) for each component of for each sequentially, starting with the coarsest resolution coefficients and progressively increasing the resolution. For each sampled value, we reconstruct the surface using Eq. (7) and (8), and evaluate the nearest neighbor energy (Eq. (12)). The parameter value that minimizes this energy is then taken as the estimate for this parameter.
Computational Complexity Since we use a sampling approach to minimize the energy for the local model, the complexity depends not on the number of iterations of a nonlinear optimization, but on the number of samples per parameter. The number of wavelet coefficients is equal to the number of vertices in the subdivision mesh, which we denote by . Note, however, that increases as the number of subdivision levels increases. In our experiments, we resample all training shapes with vertices. The complexity of the nearest neighbor search remains unchanged, as does the cost of evaluating the energy. However, the energy must be evaluated times for each of the coefficients. Hence, the overall complexity is , which is dominated by the part. Thus, assuming , we expect the local model, with its higherdimensional representation to take much longer to fit to the same point cloud. However, a tradeoff between detail and running time can be made by fitting coefficients only up to some level less than the number of levels in the wavelet decomposition.
5 Comparative Evaluation
In this section we evaluate both the fitting speed and quality of the global and local models. We evaluate both models for different initial alignment strategies and for different types of severe occlusion. Furthermore, we evaluate the models qualitatively for noisy scan data acquired using affordable stereo and range camera setups.
5.1 Experimental Setup
Training Data For training, we use the neutral expressions of subjects from the BU3DFE database [67]. This database contains relatively clean surfaces without occlusions, and a typical cropped face contains about 7500 vertices. Furthermore, each cropped face is equipped with 83 landmark points.
Parameterization We parametrize the database using the method of Salazar et al. [60] that deforms a template to each input face. This method is capable of predicting landmark points to aid in the template fitting. However, since we are given manually placed landmarks, our algorithm uses these instead of predicted ones. This removes a potential source of error during registration. The template we use contains 5996 vertices. We choose this lowresolution template for parametrization since the database has low resolution and does not contain small shape details. While the BU3DFE database contains six additional expressions in four different levels, we consider only neutral expressions for our comparison. The resulting registration is of highquality, which has been verified by manually inspecting each registered face.
Test Data We use a subset of 20 subjects (10 female and 10 male) of the Bosphorus database [68] to test our algorithm. Each subject is present in five occlusion levels: without occlusion, with glasses, with an occlusion of one eye by a hand, with an occlusion of the mouth by a hand, and with an occlusion of parts of the face by hair. Examples for each occlusion class can be seen in the left column of Figure 11. We chose this database as it allows the evaluation of different methods in the presence of severe occlusion. The resolution of this database is fairly high, and a typical face contains about 35000 vertices.
Each face is annotated with up to 22 landmarks. Figure 5a shows a model of the Bosphorus database with 22 landmark positions. The landmarks shown in red are used to compute an initial alignment of the test face to the learned shape space and the landmarks shown in green are used for error evaluation. Not all of the landmarks may be present in the database for two reasons: first, landmarks may be missing due to occlusion, and second, some landmarks are placed erroneously and we manually removed these landmarks. The landmarks that are present are not perfectly located, as shown in Figure 5b. Here, the location of the landmark at the nose tip is slightly shifted for two models of the same identity. We observed that in our test set the location of the landmarks usually varies by as much as across different scans.
Implementation Details In all experiments, we set the nearest neighbor distance truncation threshold to , and we set the parameter controlling the size of the hyperbox prior to . For the global model, we use 30 principal components. For the local model, we use a base grid of size and levels of subdivision.
Timing The global model can be fitted in under a minute per face for the data we are using. For the wavelet model the more levels that are used, the more accurate the reconstruction becomes, but also the more timeconsuming. While the reconstruction up to level zero runs in a few seconds, the reconstruction up to level five runs in slightly over one hour. Figure 7 shows the reconstruction of a Bosphorus scan when optimizing the shape coefficients of the local shape space up to different levels. This gives a way to trade off computation time and reconstruction accuracy. For all the experiments in the following, we evaluate the accurate reconstructions up to level 5.
level 0  level 1  level 2  level 3  level 4  level 5  input data 
5.2 Error Measures
In addition to the statistical measures discussed in Section 3.1 to evaluate the shape space itself, we use three ways to evaluate the fitting results for different initialization strategies and under different types of occlusion. The first two evaluation measures are quantitative measures, while the third measure is qualitative.
Statistical Measures A shape space should ideally be compact, general, and specific. Figure 8 shows that for the global PCA model, principal components explain more than of the data variability. Furthermore, for more than components, the generalization error only decreases slightly, which implies that the benefit of choosing more components is small. Finally, the specificity error still increases for more than components, which means the model represents plausible faces. Hence, we choose dimensionality for the global model. In the last column of Figure 8, random samples are chosen and the mean and standard deviation over all samples are shown.
For the Wavelet model, compactness and generalization are predetermined by the fact that we retain all variability in the model. For compactness, the percentage of variability explained by the model is . For generalization, the error will be . For specificity, we measured the mean point distance to be with a standard deviation of . As expected, this is slightly higher than the value for the global model, since the local wavelet model is less specific to the training data.
Fitting 10fold Cross Validation We perform a fold cross valiation on the training data, to evaluate the generalization of the model fitting for the global and local models. We first randomly split the training data into groups, where each group consists of half male and half female subjects. We then learn a statistical shape model on groups and fit it to each face of the remaining group. Figure 9 shows the error for both models, measured by the distance between each vertex of the fitting result and its corresponding vertex on the input face. Both models fit the data well, where of the vertices are with error for the global model and for the local model. The error for the global method is slightly smaller than for the local method.
Landmark Distance We evaluate the fitting quality using a subset of the landmarks (the ones shown in blue in Figure 5a) located on the input data. Note that these landmarks are not used in the initial alignment. Also recall that not all of these landmarks are present in all target scans. We only evaluate the error for those landmarks present in a given scan. The landmarks that are present in the test data are considered the ground truth landmark locations. We manually placed all of the landmarks used for testing on the mean shape of the aligned training data. The position of these landmarks after fitting are the estimated landmark positions. The distance between these estimates and the available ground truth landmark positions is the error we measure for the test data. This evaluation is commonly considered an accurate form of error measurement. However, unfortunately, it is possible only for a small subset of surface points since we require ground truth landmarks for this test. In the following, we give a less accurate, yet more dense evaluation, using a surface distance.
Surface Distance We evaluate the distance between the input scan and the fitted model over the entire surface by computing the distance to the nearest neighbor on the input data for each point on the fitted model. This gives a lower bound on the fitting error in terms of semantically meaningful correspondences, but it can be computed for the entire surface.
Visual Qualitative Evaluation Finally, we evaluate the results visually by showing the distance from the fitted model to the input scan.
5.3 Influence of Initialization
We first evaluate the influence of the two different initialization strategies discussed in Section 4.1 on the results. Tables 2 and 3 give the error statistics of the landmark distance and the surface distance for models without occlusions. Note that the two initialization strategies yield results of similar quality for our test database for both statistical models. This implies that both models are robust with respect to changes in the initialization.
We observed a similar behavior of the two initialization strategies for all occlusion levels with the exception that for two models with severe hair occlusion the initialization using Spin images fails, which leads to poor results.
Initialization  Mean  Median  Std. Dev.  Max 

Global Model  
manual landmarks  5.74  4.75  4.15  26.39 
Spin image alignment  5.51  4.50  3.90  25.39 
Local Model  
manual landmarks  5.96  5.09  4.29  30.07 
Spin image alignment  5.64  4.50  4.00  25.27 
Initialization  Mean  Median  Std. Dev.  Max 

Global Model  
manual landmarks  1.06  0.81  0.91  14.00 
Spin image alignment  0.91  0.74  0.64  9.18 
Local Model  
manual landmarks  0.80  0.47  1.06  15.63 
Spin image alignment  0.63  0.43  0.72  16.51 
5.4 Influence of Occlusion
We now evaluate the methods in the presence of severe occlusion. For these results, we initialize using the given landmarks to remove a potential source of error. Table 4 shows the statistics of the landmark distances. Recall that both methods are evaluated against ground truth landmarks that may be shifted by up to , as shown in Figure 5. In light of this, the quality of the reconstruction from the two models, as measured by landmark distance, is not distinguishable.
Model  Occlusion  Mean  Median  Std. Dev.  Max 

global  none  5.74  4.75  4.15  26.39 
glasses  6.68  5.64  4.71  31.05  
eye  6.87  5.29  5.29  35.74  
mouth  8.83  6.65  6.73  38.48  
hair  6.88  5.43  5.07  35.08  
local  none  5.96  5.09  4.29  30.07 
glasses  6.75  5.48  4.83  30.48  
eye  6.99  5.71  5.47  35.87  
mouth  7.30  6.01  4.83  27.62  
hair  6.79  5.34  5.28  36.56 
Table 5 shows the surface distances. We see that the local model produces slightly higher mean and standard deviation, but lower median errors. This reflects the fact that for the local model there are a relatively few points where there is not enough pull from the energy function to the data, and these points are not fitted well. In areas of the surface where the data is close enough to the initial alignment, however, the local model better fits to the surface than the global model due to the retained variability in the model. Conversely, the global model does not get as close to any localized area, but the global information allows the overall shape to guide it in areas where the initial alignment is not close enough to the data.
Model  Occlusion  Mean  Median  Std. Dev.  Max 

global  none  1.06  0.81  0.91  14.00 
glasses  1.31  0.97  1.17  13.14  
eye  2.69  1.16  5.01  65.63  
mouth  3.57  1.62  5.33  46.93  
hair  3.17  1.34  5.63  46.53  
local  none  0.80  0.47  1.06  15.63 
glasses  1.10  0.59  1.46  16.73  
eye  2.03  0.56  4.51  64.79  
mouth  2.42  0.75  4.13  47.70  
hair  2.28  0.67  4.80  51.26 
Figure 10 shows the color coded median surface distances for all points. For the results without occlusion, in most regions of the face, the local shape space yields results that are closer to the input surface. However, at the nose tip and the chin, the global model is closer to the input surface than the local model. The reason is that the initialization is often poor in these regions and that as a result the local model does not fit these areas to the surface. For the models with occlusion, the additional error caused by the occlusion is generally more localized when using the local shape space than when using the global shape space. This is especially visible in the region around the left eye for the examples where the right eye is occluded by a hand (third row of Figure 10). For the local shape space, the region around the left eye has low average fitting error, while for the global shape space, this region has larger average fitting error because it is influenced by errors in the (symmetric) region around the right eye.
no occlusion  
glasses  
eye  
mouth  
hair  
Global  Local 
Finally, we show a visual evaluation. Figure 11 shows some examples of the fitted models for visual evaluation. Both models fit the shape model close to the input data for all of the examples. Note that overall, the results of the local method capture more shape detail than the results of the global method and that in most areas of the face, the results of the local method are fitted closer to the data than the results of the global method. A notable exception is the nose area of the subject shown in the last row of the figure. The reason is that the initialization is poor in this region, which is discussed in detail above.
The third row of Figure 11 shows a facial expression that is asymmetric in the cheek area. The output of the global method is a fairly symmetric face since the global shape prior learned the symmetric structure of the face. The output of the local method correctly captures the asymmetry in the reconstruction since the local shape prior allows for more flexibility in localized shape differences.
no occlusion  

glasses  
eye  
mouth  
hair  
input  global model  local model 
Figure 12 shows the results for a challenging case where a large part of the input face is occluded by hair. Note that in spite of the large occlusions, visually satisfactory results are found by both methods.
input  global model  local model 
5.5 Evaluation with Noisy Scans
Figure 13 shows two results obtained using noisy and incomplete stereo and range data. The 3D stereo data used as input to our comparison is obtained using the approach by Brunton et al. [69] from two input images. The resulting point cloud has missing data, which is typical for data obtained using passive stereo approaches. The range data used as input to our comparison is obtained using a Kinect sensor. This dataset has low resolution, missing data, and significant data noise. In spite of these problems, both models fit the shape to the input data well in both cases. As in previous experiments, the result using the global model contains less localized shape detail than the result using the local model.
input  global model  local model 
5.6 Influence of Fitting Parameters
Values of fitting parameters used for both models were held constant in this study. This includes the threshold for nearest neighbor distance used in Eq. (12) to limit the influence of noise, occlusions, or other outliers in pulling the model to the wrong solution. Another example is the number of standard deviations model parameters are allowed to be from the mean, as controlled by in Eq. (14).
As mentioned in Section 5.1, in this study the settings and were used. These parameters influence the two models in different ways. For example, allowing greater deviation from the mean () would allow the global model to capture more detail present in the input. Since both models treat the model parameters as independent, and the global model has many fewer parameters than the local one, setting all parameters to the same distance from the mean results in a much higher probability in the global case than in the local one. On the other hand, allowing a higher threshold for nearest neighbor distance () would allow the local model to fit better to the nose in some scans. Due to the (leastsquares) rigid initial alignment, the tip and ridge of the nose of the model template are often beyond this threshold from the input point cloud, and the local influence of the model parameters means the shape of the rest of the face cannot pull this part of the model towards the data.
The parameters specific to each model were chosen based on the application to faces, or from the training data. The number of principle components kept in the global model was chosen based on the measures of the statistical model given in Section 3.1. The dimensions of the base grid of the local model was chosen so that the a single scaling coefficient was assigned to the approximate area of important facial features. Choosing a coarser base grid might allow global shape properties to be better reconstructed. The number of samples in the optimization of the local model was chosen to balance precision versus running time.
6 Implications and Observations
Let us now consider possible applications, and how this study can provide insight into which statistical model to use for a particular application. We have seen that increased localization and decorrelation in the wavelet model allow better generalization and therefore better detail recovery in the fitting process. Such detail can be important in applications such as ergonomic design (eg. designing sizes for eyeglasses that sufficiently cover a population from a relatively small training set), telepresence, and detailed facial capture in the entertainment industry, where state of the art methods require subject specific models [70] or sophisticated setups [71]. On the other hand, less detail and a lowerdimensional model are generally preferable for applications such as face and expression recognition [72, 73], modelbased manipulation of images and videos [7, 12, 13, 14], and control of virtual avatars [74].
By making our statistical models public [1], along with code to use them, we make it possible for researchers with knowledge of specific applications to test them thoroughly as they see fit, for their intended application. The provided code uses the barebones fitting energy we have used in this paper to illuminate the differences between the two types of models. However, it should be straightforward for others to modify the code to experiment with adding more terms to this energy, such as smoothing terms, landmark terms, or descriptor matching terms.
An interesting observation in Section 5.3 is that performing initial alignment with the spinimage + RANSAC approach [61] results in equal or better surface fitting and landmark errors than using manual landmarks. This means that for neutral expression, rather wellknown methods suffice to rigidly align face scans to a model. It also likely reflects the fact that landmarks in the test database [68], were placed digitally by clicking the scan rather than placed physically by an anthropometrist.
7 Conclusion
In this paper we have reviewed different statistical shape models in the focused context of 3D data fitting, and performed a comparative analysis, both theoretically and experimentally, of global and local statistical shape models for fitting to 3D face data. We have found the following differences between the two types of models: Local models capture details better at the cost of greater computational requirements. This is in part due to the optimization strategy used in this investigation (a sampling strategy that avoids local minima), but is also due to the much higher dimensionality of the local model. The global model has much lower dimensionality and can thus be fitted to input data much faster. In some cases, the global model better captures the overall shape, height and width, of the face. The local model avoids overfitting, because local surface patches are not likely to be biased for a particular database the way the shape of the entire face can be; local surface patches from human faces have much lower shape variation than entire faces, hence a limited training set has a better chance of capturing the full variability in the local model. The local model also better contains erroneous reconstruction due to occlusion to the affected areas, whereas the global model typically captures approximately symmetric shapes of human faces; an occlusion of the left side will cause poor fitting on the right as well. The local model can capture additional details by subdivision resampling.
While to a large extent, these conclusions reflect the motivations behind the use of these models, our observations about the dimensionality of the shape space with respect to the size of the training set provide useful insight into the appropriate choice of model for practitioners. If a limited number of training samples are available, the local wavelet model is preferable because each training sample will be decomposed into many independent lowdimensional samples. Conversely, if a vast amount of parametrized data is available, the global model may be preferable. If for a particular class of shapes, a partbased decomposition can be reliably and efficiently obtained, a partbased model may be preferable. This is reflected in how different models have been preferred in the literature for different types of data: for human faces and bodies, global and partbased models have traditionally been preferred, whereas for medical data, local models are more commonly used.
Applications of statistical model fitting to 3D data include face or body recognition, expression or pose recognition, biometric passwords, and virtual change rooms. In this study, we have kept the focus on model fitting to static data, however, tracking of noisy dynamic point clouds with occlusion remains an open problem.
Looking forward, sparse statistical models are a rapidly developing and expanding research area, and the use of sparsity inducing priors in statistical model fitting of 3D data provides many open avenues for future research. Initial steps have been taken in this direction [75], although so far this is limited to sparsity through locality. Recent methods for biological data formulate finding the basis of the shape space and attaching the statistical prior as an optimization problem that allows to tradeoff between locality and compactness [6]. Further, simultaneous parametrization and model learning has only been addressed by a few existing methods [42, 33]. There is no doubt room for more innovation in this direction.
Acknowledgments
We thank Eric Dubois, Jochen Lang, Chang Shu, Michael Wand, and Tino Weinkauf for helpful discussions. We thank the anonymous reviewers for their helpful and insigntful feedback and suggestions. This work has been partially funded by the Cluster of Excellence on Multimodal Computing and Interaction within the Excellence Initiative of the German Federal Government.
References

[1]
T. Bolkart, A. Brunton, A. Salazar, S. Wuhrer,
Statistical 3d
shape models of human faces (2013).
URL http://statisticalfacemodels.mmci.unisaarland.de/  [2] H. Li, T. Weise, M. Pauly, Examplebased facial rigging, ACM Transactions on Graphics 29 (4) (2010) 32:1–6.
 [3] O. van Kaick, H. Zhang, G. Hamarneh, D. CohenOr, A survey on shape correspondence, Computer Graphics Forum 30 (6) (2011) 1681–1707.
 [4] G. Tam, Z.Q. Cheng, Y.K. Lai, F. Langbein, Y. Liu, D. Marshall, R. Martin, X.F. Sun, P. Rosin, Registration of 3d point clouds and meshes: A survey from rigid to nonrigid, IEEE Transactions on Visualization and Computer Graphics 19 (7) (2013) 1199–1217.
 [5] T. Cashman, A. Fitzgibbon, What shape are dolphins? Building 3d morphable models from 2d images, IEEE Transactions on Pattern Analysis and Machine Intelligence 35 (2013) 232–244.
 [6] D. A. Alcantara, O. Carmichael, W. HarcourtSmith, K. Sterner, S. R. Frost, R. Dutton, P. Thompson, E. Delson, N. Amenta, Exploration of shape variation using localized components analysis, IEEE Transactions on Pattern Analysis and Machine Intelligence 31 (8) (2009) 1510–1516.
 [7] V. Blanz, T. Vetter, A morphable model for the synthesis of 3d faces, in: ACM Conference on Computer Graphics and Interactive Techniques, 1999, pp. 187–194.
 [8] B. Amberg, A. Blake, A. Fitzgibbon, S. Romdhani, T. Vetter, Reconstructing high quality facesurfaces using model based stereo, in: IEEE International Conference on Computer Vision, 2007, pp. 1–8.

[9]
A. Patel, W. Smith, 3d morphable face models revisited, in: IEEE Conference on Computer Vision and Pattern Recognition, 2009, pp. 1327–1334.
 [10] F. Yang, J. Wang, E. Shechtman, L. Bourdev, D. Metaxas, Expression flow for 3daware face component transfer, ACM Transactions on Graphics 30 (4) (2011) 60:1–10.
 [11] B. Amberg, R. Knothe, T. Vetter, Expression invariant 3D face recognition with a morphable model, in: IEEE International Conference on Automatic Face and Gesture Recognition, 2008, pp. 1–6.
 [12] D. Vlasic, M. Brand, H. Pfister, J. Popović, Face transfer with multilinear models, ACM Transactions on Graphics 24 (3) (2005) 426–433.
 [13] K. Dale, K. Sunkavalli, M. Johnson, D. Vlasic, W. Matusik, H. Pfister, Video face replacement, ACM Transactions on Graphics 30 (6) (2011) 130:1–10.
 [14] F. Yang, L. Bourdev, J. Wang, E. Shechtman, D. Metaxas, Facial expression editing in video using a temporallysmooth factorization, in: IEEE International Conference on Computer Vision and Pattern Recognition, 2012, pp. 861–868.
 [15] T. Bolkart, S. Wuhrer, Statistical analysis of 3d faces in motion, in: IEEE International Conference on 3D Vision, 2013, pp. 103–110.
 [16] C. Basso, A. Verri, Fitting 3d morphable models using implicit representations, in: Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications, 2007, pp. 45–52.
 [17] F. ter Haar, R. Veltkamp, 3d face model fitting for recognition, in: European Conference on Computer Vision, 2008, pp. 652–664.
 [18] M. Smet, L. V. Gool, Optimal regions for linear modelbased 3d face reconstruction, in: Asian Conference on Computer Vision, 2010, pp. 276–289.
 [19] I. Kakadiaris, G. Passalis, G. Toderici, M. Murtuza, Y. Lu, N. Karamelpatzis, T. Theoharis, Threedimensional face recognition in the presence of facial expressions: An annotated deformable model approach, IEEE Transactions on Pattern Analysis and Machine Intelligence 29 (4) (2007) 640–649.
 [20] A. Golovinskiy, W. Matusik, H. Pfister, S. Rusinkiewicz, T. Funkhouser, A statistical model for synthesis of detailed facial geometry, ACM Transactions on Graphics 25 (3) (2006) 1025–1034.
 [21] A. Brunton, C. Shu, J. Lang, E. Dubois, Wavelet modelbased stereo for fast, robust face reconstruction, in: Canadian Conference on Computer and Robot Vision, 2011, pp. 347–354.
 [22] B. Allen, B. Curless, Z. Popović, The space of human body shapes: reconstruction and parameterization from range scans, ACM Transactions on Graphics 22 (3) (2003) 587–594.
 [23] H. Seo, Y. I. Yeo, K. Wohn, 3d body reconstruction from photos based on range scan, Technologies for ELearning and Digital Entertainment (2006) 849–860.
 [24] Y. Chen, R. Cipolla, Learning shape priors for single view reconstruction, in: IEEE International Workshop on 3D Digital Imaging and Modeling, 2009, pp. 1425–1432.
 [25] J. Boisvert, C. Shu, S. Wuhrer, P. Xi, Threedimensional human shape inference from silhouettes: Reconstruction and validation, Machine Vision and Applications 24 (1) (2013) 145–157.
 [26] S. Wuhrer, C. Shu, Estimating 3d human shapes from measurements, Machine Vision and Applications 24 (6) (2013) 1133–1147.
 [27] D. Anguelov, P. Srinivasan, D. Koller, S. Thrun, J. Rodgers, J. Davis, Scape: shape completion and animation of people, ACM Transactions on Graphics 24 (3) (2005) 408–416.
 [28] P. Guan, A. Weiss, A. O. Balan, M. J. Black, Estimating human shape and pose from a single image, in: International Conference on Computer Vision, 2009, pp. 1381–1388.
 [29] A. Balan, M. Black, The naked truth: Estimating body shape under clothing, in: European Conference on Computer Vision, 2008, pp. 15–29.
 [30] S. Zhou, H. Fu, L. Liu, D. CohenOr, X. Han, Parametric reshaping of human bodies in images, ACM Transactions on Graphics 29 (4) (2010) 126:1–10.
 [31] A. Jain, T. Thormählen, H.P. Seidel, C. Theobalt, MovieReshape: Tracking and reshaping of humans in videos, ACM Transactions on Graphics 29 (5) (2010) 148:1–10.
 [32] A. Weiss, D. Hirshberg, M. Black, Home 3d body scans from noisy image and range data, in: International Conference on Computer Vision, 2011, pp. 1951–1958.
 [33] D. Hirshberg, M. Loper, E. Rachlin, M. Black, Coregistration: Simultaneous alignment and modeling of articulated 3d shape, in: European Conference on Computer Vision, 2012, pp. 242–255.
 [34] N. Hasler, C. Stoll, M. Sunkel, B. Rosenhahn, H.P. Seidel, A statistical model of human pose and body shape, Computer Graphics Forum 28 (2) (2009) 337–346.
 [35] N. Hasler, C. Stoll, B. Rosenhahn, T. Thormählen, H.P. Seidel, Estimating body shape of dressed humans, Computers & Graphics 33 (3) (2009) 211–216.
 [36] N. Hasler, H. Ackermann, B. Rosenhahn, T. Thormählen, H.P. Seidel, Multilinear pose and body shape estimation of dressed subjects from image sets, in: IEEE Conference on Computer Vision and Pattern Recognition, 2010, pp. 1823–1830.
 [37] S. Wuhrer, C. Shu, P. Xi, Postureinvariant statistical shape analysis using laplace operator, Computers & Graphics 36 (2012) 410–416.
 [38] P. Xi, W.S. Lee, C. Shu, A datadriven approach to human body cloning using a segmented body database, in: Pacific Graphics, 2007, pp. 139–147.
 [39] Y. Chen, Z. Liu, Z. Zhang, Tensorbased human body modeling, in: IEEE International Conference on Computer Vision and Pattern Recognition, 2013, pp. 105–112.
 [40] T. Cootes, C. Taylor, D. Cooper, J. Graham, Active shape models – their training and application, Computer Vision and Image Understanding 61 (1) (1995) 38–59.
 [41] T. Cootes, C. Taylor, Statistical models of appearance for medical image analysis and computer vision, in: SPIE Medical Imaging, 2001, pp. 236–248.
 [42] R. Davies, C. Twining, C. Taylor, Statistical Models of Shape: Optimisation and Evaluation, Springer, 2008.
 [43] P. Fletcher, C. Lu, S. Pizer, S. Joshi, Principal geodesic analysis for the study of nonlinear statistics of shape, IEEE Transactions on Medical Imaging 23 (8) (2004) 995–1005.
 [44] M. Toews, D. Collins, T. Arbel, A statistical partsbased appearance model of intersubject variability, in: International Conference on Medical Image Computing and Computer Assisted Intervention, 2006, pp. 232––240.
 [45] F. Lecron, J. Boisvert, S. Mahmoudi, H. Labelle, M. Benjelloun, Fast 3d spine reconstruction of postoperative patients using a multilevel statistical model, in: International Conference on Medical Image Computing and Computer Assisted Intervention, 2012, pp. 446–453.
 [46] C. Davatzikos, X. Tao, D. Shen, Hierarchical active shape models, using the wavelet transform, IEEE Transactions on Medical Imaging 22 (3) (2003) 414–423.
 [47] D. Nain, S. Haker, A. Bobick, A. Tannenbaum, Multiscale 3d shape analysis using spherical wavelets, in: International Conference on Medical Image Computing and Computer Assisted Intervention, 2005, pp. 459–467.
 [48] D. Nain, S. Haker, A. Bobick, A. Tannenbaum, Shapedriven 3d segmentation using spherical wavelets, in: International Conference on Medical Image Computing and Computer Assisted Intervention, 2006, pp. 66–74.
 [49] Y. Li, T.S. Tan, I. Volkau, W. Nowinski, Modelguided segmentation of 3D neuroradiological image using statistical surface wavelet model, in: IEEE International Conference on Computer Vision and Pattern Recognition, 2007, pp. 1–7.
 [50] S. Essafi, G. Langs, Hierarchical 3D diffusion wavelet shape priors, in: IEEE International Conference on Computer Vision, 2009, pp. 1717–1724.
 [51] P. Yu, B. T. T. Yeo, P. E. Grant, B. Fischl, P. Golland, Cortical folding development study based on overcomplete spherical wavelets, in: IEEE International Conference on Computer Vision, 2007, pp. 1–8.
 [52] T. Weise, S. Bouaziz, H. Li, M. Pauly, Realtime performancebased facial animation, ACM Transactions on Graphics 30 (4) (2011) 77:1–10.
 [53] I. Dryden, K. Mardia, Statistical Shape Analysis, Wiley, 2002.
 [54] M. Styner, K. Rajamani, L.P. Nolte, G. Zsemlye, G. Szekely, C. Taylor, R. Davies, Evaluation of 3d correspondence methods for model building, in: Information Processing in Medical Imaging, 2003, pp. 63–75.
 [55] W. Sweldens, The lifting scheme: A customdesign construction of biorthogonal wavelets, Applied and Computational Harmonic Analysis 3 (2) (1996) 186–200.
 [56] S. Mallat, A Wavelet Tour of Signal Processing, Elsevier, 1999.
 [57] P. Schröder, W. Sweldens, Spherical wavelets: Efficiently representing functions on the sphere, in: ACM Conference on Computer Graphics and Interactive Techniques, 1995, pp. 161–172.
 [58] M. Bertram, M. Duchaineau, B. Hamann, K. I. Joy, Generalized BSpline subdivisionsurface wavelets for geometry compression, IEEE Transactions on Visualization and Graphics 10 (3) (2004) 326–338.
 [59] C. Creusot, N. Pears, J. Austin, 3D face landmark labelling, in: ACM Workshop on 3D Object Retrieval, 2010, pp. 27–32.
 [60] A. Salazar, S. Wuhrer, S. Chu, F. Prieto, Fully automatic expressioninvariant face correspondence, Machine Vision and Applications 25 (4) (2014) 859–879.
 [61] A. Johnson, M. Hebert, Recognizing objects by matching oriented points, in: Conference on Computer Vision and Pattern Recognition, 1997, pp. 684–692.
 [62] M. Fischler, R. Bolles, Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography, Communications of the ACM 24 (6) (1981) 381–395.

[63]
D. Mount, S. Arya, ANN: A library
for approximate nearest neighbor searching (2010).
URL http://www.cs.umd.edu/~mount/ANN/  [64] A. Patel, W. Smith, Exploring the identity manifold: Constrained operations in face space, in: European Conference on Computer Vision, 2010, pp. 112–125.
 [65] D. Liu, J. Nocedal, On the limited memory method for large scale optimization, Mathematical Programming 45 (3) (1989) 503–528.
 [66] D. Lee, C. Wong, Worstcase analysis for region and partial region searches in multidimensional binary search trees and balanced quad trees, Acta Informatica 9 (1) (1977) 23–29.
 [67] L. Yin, X. We, Y. Sun, J. Wang, M. Rosato, A 3D facial expression database for facial behavior research, in: IEEE International Conference on Automatic Face and Gesture Recognition, 2006, pp. 211–216.
 [68] A. Savran, N. Alyuz, H. Dibeklioglu, O. Celiktutan, B. Gökberk, B. Sankur, L. Akarun, Bosphorus database for 3D face analysis, in: Workshop on Biometrics and Identity Management, 2008, pp. 47–56.
 [69] A. Brunton, J. Lang, E. Dubois, Efficient multiscale stereo of highresolution planar and spherical images, in: IEEE Conference on 3D Imaging, Modeling, Processing, Visualization, and Transmission, 2012, pp. 120–127.
 [70] P. Garrido, L. Valgaerts, C. Wu, C. Theobalt, Reconstructing detailed dynamic face geometry from monocular video, ACM Transactions on Graphics 32 (6) (2013) 158:1–10.
 [71] T. Beeler, F. Hahn, D. Bradley, B. Bickel, P. Beardsley, C. Gotsman, R. W. Sumner, M. Gross, Highquality passive facial performance capture using anchor frames, ACM Transactions on Graphics 30 (2011) 75:1–10.
 [72] K. Scherbaum, J. Petterson, R. S. Feris, V. Blanz, H.P. Seidel, Fast face detector training using tailored views, in: IEEE International Conference on Computer Vision, 2013.
 [73] I. Mpiperis, S. Malassiotis, M. G. Strintzis, Bilinear models for 3d face and facial expression recognition, IEEE Transactions on Information Forensics and Security 3 (3) (2008) 498–511.
 [74] H. Li, J. Yu, Y. Ye, C. Bregler, Realtime facial animation with onthefly correctives, ACM Transactions on Graphics 32 (4).
 [75] T. Neumann, K. Varanasi, S. Wenger, M. Wacker, M. Magnor, C. Theobalt, Sparse localized deformation components, ACM Transactions on Graphics 32 (6) (2013) 179:1–10.