Data-Driven Physical Face Inversion

TL;DR

This paper introduces a simple finite element simulation method for facial animation, along with a novel technique to automatically recover actor-specific physical parameters from minimal 3D head pose data, enabling realistic and personalized face simulations.

Contribution

It presents a new approach to recover physical face parameters using minimal data, improving the realism and ease of tuning in physical face animation methods.

Findings

The method accurately matches real actor face motion.

Recovered parameters enable realistic, actor-specific face simulations.

Analysis of spatial material layouts improves simulation quality.

Abstract

Facial animation is one of the most challenging problems in computer graphics, and it is often solved using linear heuristics like blend-shape rigging. More expressive approaches like physical simulation have emerged, but these methods are very difficult to tune, especially when simulating a real actor's face. We propose to use a simple finite element simulation approach for face animation, and present a novel method for recovering the required simulation parameters in order to best match a real actor's face motion. Our method involves reconstructing a very small number of head poses of the actor in 3D, where the head poses span different configurations of force directions due to gravity. Our algorithm can then automatically recover both the gravity-free rest shape of the face as well as the spatially-varying physical material stiffness such that a forward simulation will match the captured targets as closely as possible. As a result, our system can produce actor-specific, physical parameters that can be immediately used in recent physical simulation methods for faces. Furthermore, as the simulation results depend heavily on the chosen spatial layout of material clusters, we analyze and compare different spatial layouts.

Figures

Tables

Table 1. Table 1: Validation results. The validation sequence consists of 136 poses, out of which three are used as optimization targets. The error is measured as the total absolute point to point vertex error between the simulation result and the tracked mesh.

Configuration	Total Error (mm)	Mean Vertex Error (mm)
Naive	263048	1.31
Layout 1	210186	1.05
Layout 2	212624	1.06

Equations

$$W(\mathbf{X},\mathbf{x},\mathbf{P})=\sum _e{W}_e(\mathbf{X},\mathbf{x},\mathbf{P}),$$

$${\nabla }_{\mathbf{x}}W(\mathbf{X},\mathbf{x},\mathbf{P})={\mathbf{f}}_{\text{ext}},$$

$$g(\mathbf{X},\mathbf{P})=\frac{1}{2}{\Vert \mathbf{Sx}(\mathbf{X},\mathbf{P})-{\bar{\mathbf{x}}}_s\Vert }^2+\alpha R(\mathbf{X}),$$

$$\text{subject to}{\mathbf{P}}_l\le \mathbf{P}\le {\mathbf{P}}_u,$$

$$\left[\begin{array}[]{c}\frac{\text{d}g(\mathbf{X},\mathbf{P})}{\text{d}\mathbf{X}}\\ \frac{\text{d}g(\mathbf{X},\mathbf{P})}{\text{d}\mathbf{P}}\end{array}\right]=\left[\begin{array}[]{c}\alpha\frac{dR(\mathbf{X})}{d\mathbf{X}}+\left(\mathbf{Sx}(\mathbf{X,P})-\bar{\mathbf{x}}_{s}\right)^{T}\mathbf{S}\frac{\text{d}\,\mathbf{x}(\mathbf{X},\mathbf{P})}{\text{d}\,\mathbf{X}}\\ \left(\mathbf{Sx}(\mathbf{X,P})-\bar{\mathbf{x}}_{s}\right)^{T}\mathbf{S}\frac{\text{d}\,\mathbf{x}(\mathbf{X},\mathbf{P})}{\text{d}\,\mathbf{P}}\end{array}\right].$$

$$-\frac{\partial {\mathbf{f}}_{\text{ext}}(\mathbf{y})}{\partial \mathbf{y}}+\frac{{\partial }^{2}W(\mathbf{x},\mathbf{y})}{\partial \mathbf{y}\partial \mathbf{x}}+\frac{{\partial }^{2}W(\mathbf{x},\mathbf{y})}{\partial {\mathbf{x}}^2}\frac{\text{d}\mathbf{x}(\mathbf{y})}{\text{d}\mathbf{y}}\stackrel{!}{=}\mathbf{O},$$

$$\frac{{\partial }^{2}W(\mathbf{x},\mathbf{y})}{\partial {\mathbf{x}}^2}\frac{\text{d}\mathbf{x}(\mathbf{y})}{\text{d}\mathbf{y}}=\mathbf{b}(\mathbf{x},\mathbf{y})$$

$$\mathbf{b}(\mathbf{x},\mathbf{y})=\frac{\partial {\mathbf{f}}_{\text{ext}}(\mathbf{y})}{\partial \mathbf{y}}-\frac{{\partial }^{2}W(\mathbf{x},\mathbf{y})}{\partial \mathbf{y}\partial \mathbf{x}}$$

$$\frac{\text{d}\mathbf{x}(\mathbf{y})}{\text{d}\mathbf{y}}={\left(\frac{{\partial }^{2}W(\mathbf{x},\mathbf{y})}{\partial {\mathbf{x}}^2}\right)}^{-1}\mathbf{b}(\mathbf{x},\mathbf{y})$$

$$\frac{\text{d}g(\mathbf{y})}{\text{d}\mathbf{y}}=\frac{dR(\mathbf{y})}{d\mathbf{y}}+{\left(\mathbf{Sx}-{\bar{\mathbf{x}}}_s\right)}^T\mathbf{S}{\left(\frac{{\partial }^{2}W(\mathbf{x},\mathbf{y})}{\partial {\mathbf{x}}^2}\right)}^{-1}\mathbf{b}(\mathbf{x},\mathbf{y}).$$

$${\mathbf{\lambda }}^T={\left(\mathbf{Sx}-{\bar{\mathbf{x}}}_s\right)}^T\mathbf{S}{\left(\frac{{\partial }^{2}W(\mathbf{x},\mathbf{y})}{\partial {\mathbf{x}}^2}\right)}^{-1},$$

$$\left(\frac{{\partial }^{2}W(\mathbf{x},\mathbf{y})}{\partial {\mathbf{x}}^2}\right)\mathbf{\lambda }={\mathbf{S}}^T\left(\mathbf{Sx}-{\bar{\mathbf{x}}}_s\right)$$

$$\frac{\text{d}g(\mathbf{y})}{\text{d}\mathbf{y}}=\alpha \frac{dR(\mathbf{y})}{d\mathbf{y}}+{\mathbf{\lambda }}^T\mathbf{b}(\mathbf{x},\mathbf{y}).$$

$$E_e=\sum _{i=1}^c{w}_e^i{E}_c,$$

$$\frac{dg}{d{E}_c^i}=\sum _{e=1}^m\frac{dg}{d{E}_e}{w}_e^i.$$

$$R(\mathbf{X})=W({\mathbf{X}}^0,{\mathbf{X},\mathbf{P}}^0),$$

$$g(X^{(k)}-\beta \Delta \mathbf{X})\le g(X^{(k)})-\gamma \beta \Delta {\mathbf{X}}^T\nabla g(X^{(k)}),$$

$$g(\mathbf{X},\mathbf{P})$$

Full text

††volume: 38††issue: 2

1 Introduction

Creating realistic physical effects for human faces using simulation methods in computer graphics is an extremely challenging problem. On the one hand, linear surface rigs, such as blendshapes, offer intuitive control but are not suited for physical simulation and physical effects have to be added manually. On the other hand, anatomically accurate rigs are hard to create and control, and require high expertise in both human anatomy and CG animation. Recently, simple volumetric rigs have been proposed as middleground [IKNDP16, KBB*∗*17, IKKP17], which maintain the intuitive control of blendshapes but are also suited for physical simulation. Finding suitable material properties for such rigs is, however, a challenging task, which so far is only possible in a manual trial and error process. And since physical properties change significantly between humans, due to tissue distribution, age and body mass index (BMI), they have to be re-specified anew for every character. In this paper we investigate how to directly estimate the physical properties of human faces from captured data.

One possible approach to this problem would be to measure tissue properties in-vitro for the different tissues that constitute a face, and then distribute these properties assuming some knowledge of the internal tissue distributions. Of course, usually such knowledge is not available, and even in a case where the internal structure of the tissue can be observed, i.e. from MRI or CT data, and when the anatomical structures are known, recovering material parameters for simulation (e.g. stiffness) is still challenging since in-vitro uni-axial strain-stress response measured from thin samples do not directly map to the behaviour of larger tissue structures in-vivo.

We feel a more promising approach is to measure the material properties of human soft tissue in-vivo. While there have been approaches that measure human tissue at a few sparse locations using a force probe [BBO*∗*09, PRWH*∗*18], we propose a much more convenient approach that yields dense measurements across the entire face. We capture a subject under a small set of different head poses, which naturally cause the skin tissue to deform due to the known gravity force. From these measurements we can invert both the rest-shape, i.e. the shape of the face in the absence of gravity, as well as material properties, i.e. spatially-varying material stiffness. Combined, these allow to physically simulate the face, yielding results as close as possible to the captured shapes when applying the respective gravitational forces, and extrapolating to physically plausible shapes for other forces. We show that it is essential to represent the face with spatially varying material in order to better approximate the underlying anatomy, and we evaluate several face partition layouts that were presented in the literature. Since the proposed method is agnostic to the layout, we can determine optimal simulation properties for any layout.

To summarize, we present a novel system to estimate rest-shape and spatially varying material properties of real human faces from a few captured 3D scans, making the acquisition and creation of simulation ready facial rigs very practical. While the approach is formulated in the context of human faces, the underlying theory and algorithms are generic and could find application in other areas as well.

2 Related Work

Arguably, the most common approach to animate faces is blendshape based facial animation [LAR*∗*14]. However, it is well-known that blendshapes are quite limiting in that they provide only a linear subspace. Recently, researchers have turned to physical simulation in an attempt to achieve more expressive, non-linear facial animation.

2.1 Physical Faces

One of the first methods for physics-based facial animation was presented by Sifakis et al. [SNF05], who built a detailed face rig consisting of a complete, anatomically accurate muscle structure, created manually from the actor’s medical data. Building the muscle structure for an actor is a time consuming process, therefore Cong et al. [CBJ*∗*15] developed an automatic way to transfer a template anatomy to target input faces. These transferred muscle-based rigs can then be artistically refined by modifying the tracks that muscles follow during activation [CBF16]. Ichim et al. \shortciteichim2017phace also fit a template model of bones, muscles and flesh to facial scans. Their method succeeds by solving for the muscle activation parameters that best fit the input scans during forward simulation, and thus produces an actor-specific physical face mesh for animation. A similar approach for full bodies was explored by Kadleček et al. \shortcitekadlecek2016.

Alternative to transferring template anatomy, another approach to avoid the costly setup time of anatomically accurate simulations is to use other representations of the face, for example surface-based [YSZ09, MWF*∗*12, BSC16, BS18] or volumetric rigs [LXB17, IKNDP16, KBB*∗*17]. On the side of surface representations, Ma et al. \shortcitema2012 use a mass-spring system to build a blend shape model that incorporates physical interaction. You et al. \shortciteyou2009, as well as Barrielle et al. \shortcitebarrielle2016 propose to define forces at the vertices of a face mesh and blend different forces that correspond to different face shapes in order to generate facial animation in a simulation setting, which can even be extended to real-time animation [BS18]. While surface based simulation methods are typically easier to setup, they cannot create volume based dynamics such as of soft tissue, which is one of the benefits of volumetric rigs. On the volumetric side, Li et al \shortciteli2017tvcg propose a method to enrich triangle mesh animations by fitting a tetrahedral mesh, applying physics, and then transferring the secondary motion and collision resolution to the original mesh. Ichim et al. \shortciteIchim2016, as well as Kozlov et al. \shortcitekozlov2017 propose to more intricately couple blendshape-based facial animation with user-specific volumetric rigs, creating blend-volumes that can be used in a finite-element simulation approach. The work of Ichim et al. \shortciteIchim2016 allows for several interesting physical effects but does not account for expression-specific dynamics. Kozlov et al.  \shortcitekozlov2017 focus on the creation of expression-specific physical effects, but the drawback is that spatially-varying material parameters need to be painted and set manually by an artist for each expression. Our method is complementary to Kozlov et al. \shortcitekozlov2017 as we aim to automatically determine the material parameters for a similar volumetric simulation rig.

2.2 Data-Driven Parameter Estimation

Our work falls into the category of estimating physical parameters from real world capture data. Several methods have been presented, outside the focus area of facial animation, for example on the topic of capturing cloth simulation parameters [BTH*∗*03, WOR11, MBT*∗*12], optimizing rest shape geometry of elastic objects [CZXZ14], estimating damping properties for dynamics [XB17, CLMK17], or data-driven modeling of generic soft objects [WWY*∗*15, MMO16].

Our work estimates both nonlinear heterogeneous material distribution and rest shape geometry for human faces, under a unified optimization framework of sensitivity analysis. Similar to recent methods [XLCB15, WWY*∗*15], we estimate the heterogeneous distribution of material parameters in a reduced parametric space. However, we parameterize the spatial distribution into several material clusters based on the physical facial structures whereas previous work project the distribution into smooth basis functions without correspondence to real-world data.

Notably, some approaches in the category of data-driven soft object modeling have been designed to measure human tissue. Bickel et al. \shortcitebickel2009capture fit a stress-strain curve to elastostatic deformation of skin tissue under known forces. More recently, Pai et al. \shortcitepai2018human measure the mechanical properties of the human body with a new handheld device. Kim et al. \shortcitekim2017 learn an active and passive tissue segmentation and the material parameters for the active tissue for whole-body deformations. Finally, it is worth mentioning that Pons-Moll et al. \shortcitepons-moll2015 took a non-physics based approach, but rather a learning approach to model the deformation caused by soft-tissue dynamics in an application of full body animation, by scanning over 40,000 poses of real people. A similar approach could be used for more expressive facial animation without physical simulation.

3 Overview

As illustrated in Fig. Data-Driven Physical Face Inversion, the input to our approach is a sparse set of aligned facial surface deformations under gravity, captured at different orientations of the head with a fixed expression (Section 6). Our approach outputs the simulation mesh geometry at rest (without gravity) and spatially-varying material stiffness across the deformable volume of the face, such that when we forward simulate the face model under gravity, it will match the target surface deformations.

We model the deformable face with a tetrahedral mesh and physically simulate the deformation using the Finite Element Method (FEM) with a nonlinear material model (Section 4). To obtain the physical parameters for forward simulation, we formulate an inverse optimization to match the captured facial deformations. We optimize the unknown physical parameters, by alternating between optimizing the rest shape geometry and material stiffness distribution using sensitivity analysis (Section 5). We regularize our inverse problem by aggregating the material distribution into spatial clusters (Section 5.2) and asking for minimal rest-shape geometry changes to our initial guess (Section 5.3).

In Section 7, we analyze the optimized results from different spatial material layouts. We validate our optimized physical parameters by simulating to different poses from the capture sequence, which were not part of the optimization targets.

4 Physically-Based Face Simulation

Inspired by [KBB*∗*17] we do not explicitly model the underlying anatomy but instead treat the entire face as a single volume, with spatially varying material properties that abstract the underlying anatomical complexity. This results in much simpler creation and simulation of a human face, while still achieving physically plausible effects as recently demonstrated in various papers [IKNDP16, KBB*∗*17, IKKP17].

To simulate the deforming facial tissue, we use the Finite Element Method (FEM) simulation of nonlinear elastic materials. The deformable volume is modeled as a tetrahedral mesh, subject to fixed boundary conditions at the skull and jaw (Fig. 1).

The simulation model is governed by the underlying hyperelastic material law where the elastic potential energy $W$ is a summation of elemental contribution $W_{e}(\mathbf{X,x,P})$ as

[TABLE]

where $\mathbf{P}$ are the per-element material parameters and $\mathbf{X,x}\in{\mathcal{R}}^{3n}$ are the undeformed and deformed nodal positions, respectively, with $n$ denoting the number of nodes. We use a Neo-Hookean material [BW08] to model the nonlinearity in facial deformation, in which case $\mathbf{P}\in{\mathcal{R}}^{2m}$ is a collection of per-element Young’s modulus $E_{e}$ and Poisson’s ratio $\nu_{e}$, with $m$ denoting the number of tetrahedral elements. To improve the simulation robustness to element inversions, we employ invertible FEM simulation using thresholded deformation gradient [ITF04].

For quasi-static deformation, the deformed state $\mathbf{x}$ can be determined by balancing the elastic force and external forces as

[TABLE]

where $\mathbf{f}_{\text{ext}}\in{\mathcal{R}}^{3n}$ can be gravity force, contact or any other external loads. In section 5, we will use the quasi-static simulation for our inverse problem, where the face is statically deformed under gravity.

5 Static Inversion

Given the face model described in Section 4, under gravity the facial tissue would be statically deformed. In this section, we describe how we invert this effect by estimating simulation-ready rest-shape geometry $\mathbf{X}$ and per-element material parameters $\mathbf{P}$ based on observations of the static facial deformation under known gravity forces. In the following section, we outline a capture procedure and processing pipeline to generate such observations (Section 6).

5.1 Inverse Problem

Given a deformed configuration and the initial simulation mesh topology, we aim at solving for the undeformed vertices $\mathbf{X}$ and material parameters $\mathbf{P}$. We formulate our optimization objective as

[TABLE]

where $\bar{\mathbf{x}}_{s}\in\mathcal{R}^{3s}$ denote the subset of deformed and observed vertices on the surface and $\mathbf{S}\in\mathcal{R}^{3s\times 3n}$ is the selection matrix to the corresponding observed vertices. $R(\mathbf{X})$ is a regularization term we will use in our rest shape geometry optimization (section 5.3). $\mathbf{P}_{l}$ and $\mathbf{P}_{u}$ are the lower and upper bounds for our material parameters respectively. Note that the relationship $\mathbf{x}(\mathbf{X,P})$ is non-trivial. Given a new rest configuration or an updated set of material parameters, we first solve a quasi-static forward simulation (Equation 2) to get the deformed configuration $\mathbf{x}$ at static equilibrium before we compare the result to the observed vertices $\bar{\mathbf{x}}_{s}$.

To minimize this objective with a quasi-Newton method with BFGS, an analytical gradient is desirable,

[TABLE]

However, there is no direct analytical expression for $\frac{\text{d}\,\mathbf{x}(\mathbf{X},\mathbf{P})}{\text{d}\,\mathbf{X}}$ or $\frac{\text{d}\,\mathbf{x}(\mathbf{X},\mathbf{P})}{\text{d}\,\mathbf{P}}$ due to the aforementioned implicit relationship $\mathbf{x}(\mathbf{X},\mathbf{P})$. Therefore we propose to use sensitivity analysis to compute the analytical gradient.

Here for the simplicity of notation, we assemble all the optimization variables into vector $\mathbf{y}=[\mathbf{X},\mathbf{P}]$. Starting from Equation 2, we take the derivative with respect to the optimization variables $\mathbf{y}$ on both sides, we then get

[TABLE]

with $\mathbf{O}$ denoting the zero matrix. We note that in our case, $\mathbf{f}_{\text{ext}}$ is the gravity force only.

We could compute $\frac{\text{d}\mathbf{x}(\mathbf{y})}{\text{d}\mathbf{y}}$ by solving the equation systems

[TABLE]

with right-hand side

[TABLE]

where $\frac{\partial^{2}W(\mathbf{x,y})}{\partial\mathbf{x}^{2}}$ is the stiffness matrix we use when solving the forward problem with Newton’s method.

Adjoint Evaluation. The dimension of the right hand side (Equation 8) is $3n\times p$, where $p$ is the number of the optimization variables. Evaluating the system 7 directly requires $p$ linear system solves. Even though we can prefactorize the left hand side matrix, this is still computationally expensive when we have a high-dimensional optimization problem. Here instead we use the adjoint method [CLPS03] to reduce this into a single system solve.

We plug the explicit expression

[TABLE]

into the expression for the gradient,

[TABLE]

By introducing a column vector $\boldsymbol{\lambda}$

[TABLE]

we can avoid solving several linear equation systems but solve the adjoint system

[TABLE]

instead, resulting in the gradient as

[TABLE]

Instead of optimizing all the variables at the same time, we will use block coordinate descent for our inverse problem, alternating between material optimization 5.2 and rest shape optimization 5.3, until convergence.

5.2 Cluster-based Material Optimization

For material optimization, we keep the rest shape geometry $\mathbf{X}$ constant and optimize per-element material parameters $\mathbf{P}$. As discussed in section 4, we use Neo-Hookean material for our forward simulation where the material stiffness is parameterized with Young’s modulus $E$ and volume preservation is controlled with Poisson’s ratio $\nu.$ However, we note that for different facial tissues, they mostly differ in the Young’s modulus but are close in Poisson’s ratio with high incompressibility. We therefore only optimize heterogeneous distribution of Young’s modulus for our inverse problem but use a constant and homogeneous distribution of Poisson’s ratio (we use $\nu=0.43$).

Optimizing element-wise Young’s modulus leads to an inverse problem of $m$-dimension, where $m$ is the number of tetrahedral elements and can be very large as the simulation mesh becomes complex. In practice, high-dimensional optimization often suffers from slow converge, poor local minima and parameters overfitting. Inspired from reduced material optimization [XLCB15, WWY*∗*15] which use basis functions of smooth material distribution, here we project our element-wise material inversion problem into a spatial cluster-based optimization.

Based on anatomical structures, we first create $c$ spatial clusters ($c\ll m$) of material distribution (Section 6.5). The Young’s modulus of each tet is then a linear combination of material cluster values,

[TABLE]

where $w^{i}_{e}$ is the weight to each cluster for element $e$ and $\sum_{i=1}^{c}w^{i}_{e}=1.$ In our case, we assign each interior element exclusively to a single cluster but for elements at the boundary of clusters, we do linear interpolation between the neighboring clusters. By optimizing the material parameters of $c$ clusters only, we solve a much lower-dimensional problem where the analytical gradient to the $i$-th cluster parameter $E^{i}_{c}$ is calculated as

[TABLE]

The initial guess to our material optimization is a homogeneous distribution of Young’s modulus where $E_{c}=1.0$ MPa for all the clusters. We bound the stiffness to be within the physical range of Young’s modulus by setting $\mathbf{P}_{l}=0.001$ MPa and $\mathbf{P}_{u}=1.0$ MPa. For our constrained minimization, we use the standard interior point method and BFGS [NW06] as implemented in the KNITRO package.

5.3 Rest Shape Geometry Optimization

Although one of our observations will be a deformed surface geometry under gravity at a neutral expression (refer to Section 6), the undeformed positions of the simulation mesh $\mathbf{X}$ still remain unknown. To solve for the rest shape geometry, we keep the material parameters unchanged and optimize the rest shape $\mathbf{X}$. Note that we need to find undeformed positions for all the nodes including the interior ones, even though only a subset of the surface mesh nodes $\bar{\mathbf{x}}_{s}$ are observed by our capture system.

Similar to element-wise material optimization, the geometry optimization is also a high-dimensional problem with $3n$ number of variables. However, here we adopt a different strategy to material optimization by regularizing the objective with

[TABLE]

where $\mathbf{X}^{0}$ is the initial guess to the undeformed geometry and $\mathbf{P}^{0}$ are the some constant material parameters. The regularization term $R(\mathbf{X})$ helps maintain a consistent simulation topology by penalizing large deviation of the solution to the initial guess with the elastic energy $W$.

The initial guess to our rest shape geometry $\mathbf{X}^{0}$ is obtained by forward quasi-static simulation of the volumetric mesh by applying the opposite gravity force. We note that our initial guess is reasonably close to our solution, we therefore use gradient descent method with line search for our optimization. With the analytic gradient $\Delta\mathbf{X}=\frac{dg(\mathbf{X})}{d\mathbf{X}}$ evaluated from adjoint evaluation, we do a line search along the gradient direction. To avoid simulation mesh topology change from element inversion, we back track the maximum step length $\beta$ such that no element is inverted [SS15] and the first Wolfe condition [NW06] is satisfied:

[TABLE]

where $\gamma$ is a control parameter (we use $\gamma=10^{-4}$).

5.4 Multiple Poses

In order to robustly estimate the physical properties, we propose to include multiple poses of the same expression under different head rotations and observe the facial deformations under gravity. To compute the rest configuration and material parameters from several observations under varying pose orientations, we extend our objective as

[TABLE]

where each $o$ is an observation of the same expression in a different orientation relative to the gravity. We note that here instead of transforming the rest shape $\mathbf{X}$ according to the orientations, we align our observed targets to a single pose and rotate the direction of gravity based on the relative orientation, and perform the optimization in the canonical frame. In practice, we also observe better convergence by starting our optimization from a single pose and using the optimized result as the initial guess to our multiple-pose optimization.

6 Data Acquisition and Preparation

We now describe a procedure to capture the observations required for inversion, as described in Section 5, for a real human face.

6.1 Data Capture

Our capture setup consists of four synchronized stereo pairs of machine-vision cameras, capturing at approx. 30fps (see Fig. 2). The method of Beeler et al. \shortcitebeeler2011 is used to reconstruct geometry from performance sequences of facial motion, and generate a set of meshes in dense vertex correspondence. The face meshes are stabilized to remove the rigid motion of the head using an approximated skull position, as proposed by Beeler and Bradley \shortcitebeeler2014rigid.

6.2 Actor Performance

The capture is recorded in a standing position. The actor performs a slow rotation of his head by moving his body and keeping the head position constant relative to their neck. To prevent jaw motion relative to the skull, the actor lightly holds a small teeth guard in their mouth. The slow head rotation prevents any inertia-induced deformation, therefore the observed face deformation is assumed to be quasi-static and only due to the effects of gravity. A subset of the captured and reconstructed data is shown in Fig. 3.

6.3 Simulation Model Preparation

Laplacian deformation is used to deform a closed-surface, generic template face model to one of the target shapes. The underlying skull and jaw bones are estimated using the technique described by Zoss et al. \shortcitezoss2018.

We use Tetgen[Si06] to generate a surface conforming tetrahedralization of the volume between the face surface and the bone meshes. The correspondences between the resulting volumetric model vertices and the surface meshes of the bones are used to generate the simulation boundary constraints. The final simulation model consists of 50k vertices and 300k tetrahedral elements.

6.4 Target Preparation and Alignment

From the sequence of tracked meshes, we choose a sparse set of poses to be used as optimization targets (refer to Fig. 3). The selected targets should exhibit a wide span of poses. The effects of the soft tissue deformation exhibited by the face can be easily dominated by any alignment errors, therefore the optimization requires careful treatment of the target poses. In cases where the automatic stabilization of the meshes is visually insufficient, we manually align the head poses. The aligned frames are then used to compute a pose-specific adjusted gravity direction, as the optimization is performed in a canonical reference frame. The deformation on the vertices then serves as the observation $\mathbf{\bar{x}}_{s}$ for our optimization. See Fig. 4 for visualization of the aligned targets and computed gravity direction. Per target we generate positional constraints for about 1500 vertices of the simulation mesh, selected using a manually specified mask confined to the center of the face.

6.5 Material Cluster Preparation

The segmentation of the mesh into material clusters is done by painting colors on the surface of the volumetric mesh and propagating these colors inwards through the tetrahedral mesh’s connectivity. We evaluate two different material cluster layouts; Layout 1 is based on the material maps from Kozlov et al.\shortcitekozlov2017. This layout is comprised of 7 clusters and inherently enforces symmetry since left and right cheeks belong to the same cluster. Layout 2 is based on the one Tena et al.\shortcitetena2011interactive suggested, but we added clusters for the lips as we anticipate those to differ from the surrounding tissue, yielding a total of 19 clusters. The layouts are shown in Fig. 6.

7 Results

We validate the results quantitatively by forward simulating to a set of frames where we know the gravity vector and have ground truth surface measurements from the performance capture system. We first simulate using a naïve configuration, where a neutral scan is used as rest-shape and we manually determine a uniform material that is soft enough to produce secondary motion and other dynamic effects. While this might seem like a trivial setup, it is probably the most common approach for artists since neither geometry nor material inversion are readily available. Since this treats the entire face as a single material, results are a compromise and, for example, the nose is too soft, leading to large undesired deformations (Fig. 5, left column). Assigning spatially varying materials will remedy this problem, but picking those manually is a challenge in its own. Instead, the proposed system determines the optimal rest-shape in combination with optimal material properties for a chosen layout, which yield a lower error when used in simulation (Fig. 5, center and right columns). The overall error averaged over a total of 136 poses is summarized in Table 1.

Fig. 6 shows the two layouts used in this publication, as well as the resulting Young’s Modulus automatically determined by the proposed system. As can be seen, both the layouts yield similar semantic material assignments, i.e. the cheeks come out soft, where nose and lips are stiffer. Starting from an overall stiff material, the optimization successively updates the estimated stiffness per region until convergence as shown in Fig. 7, which is achieved after about 30 iterations. The optimized Young’s Modulus stay in the physical range for facial tissues, as reported in literature [CNJO96]. In an alternating process, we again optimize for the optimal rest shape given the resulting material configuration, which corresponds to the face under zero gravity (Fig. 8).

Performance

The simulation mesh for our face model has $50k$ vertices and $300k$ tetrahedral elements, with $7-19$ material clusters, depending on the layout. An unoptimized implementation using five optimization target poses requires up to two minutes for a single evaluation of the function on an Intel i7 3.2Ghz with eight cores, yielding a total runtime of approximately 5 hours. The energies and gradients for the different targets are evaluated in parallel.

8 Conclusion

We present a novel method to automatically estimate the physical properties required to simulate real human faces given a small number of aligned facial scans. As physical simulation for facial animation is gaining popularity, the challenge of tuning the simulation to match a real human face is becoming an increasing roadblock for using simulation approaches in several real-world scenarios (such as VFX and VR). This work represents a large step in the direction of applying simulation methods on real faces. By requiring only a sparse number of input poses, we show that the capture requirements are minimal, and we further study the effects of different spatial layouts for material clustering on the face.

Our method estimates heterogeneous material properties and rest-shape simulation geometry from multiple poses in a neutral expression. In the future, we would like to extend our method to multiple expressions, accounting for material and geometry changes across different expressions [KBB*∗*17]. Starting from our work, it is also an interesting but challenging direction to obtain physical dynamical properties from sequences of facial performance capture. For our rest shape geometry optimization we resize the step length globally to avoid element inversion, which results in robust yet slow convergence. Employing more advanced adaptive approaches [ZBK18] could lead to faster convergence and speed up the optimization. Furthermore, just like any other data-driven approach, input data noise, such as reconstruction error and geometric alignment, has an adverse effect on the optimization quality. This is aggravated by the fact that the observed deformation due to gravity is rather small, making careful geometric alignment, mouth guard, and proper instructions to maintain the same expression extremely important. In the future we would like to investigate ways to apply stronger forces to observe stronger deformations, while increasing the constraints on the acquisition procedure to reduce input noise further, leading to a better overall signal-to-noise ratio. One avenue along these lines would be to investigate sparse feature-based tracking objectives, compared to our dense point-based optimization targets. Lastly, our spatial material layout is just a coarse approximation of the facial anatomical structure and we would like to automatically devise more anatomically accurate material clustering from MRI and CT volume data.