Nonlinear Deformation Synthesis via Sparse Principal Geodesic Analysis

This paper introduces the construction of a low‐dimensional nonlinear space capturing the variability of a non‐rigid shape from a data set of example poses. The core of the approach is a Sparse Principal Geodesic Analysis (SPGA) on the Riemannian manifold of discrete shells, in which a pose of a non‐rigid shape is a point. The SPGA is invariant to rigid body motions of the poses and supports large deformation. Since the Riemannian metric measures the membrane and bending distortions of the shells, the sparsity term forces the modes to describe largely decoupled and localized deformations. This property facilitates the analysis of articulated shapes. The modes often represent characteristic articulations of the shape and usually come with a decomposing of the spanned subspace into low‐dimensional widely decoupled subspaces. For example, for human models, one expects distinct, localized modes for the bending of elbow or knee whereas some more modes are required to represent shoulder articulation. The decoupling property can be used to construct useful starting points for the computation of the nonlinear deformations via a superposition of shape submanifolds resulting from the decoupling. In a preprocessing stage, samples of the individual subspaces are computed, and, in an online phase, these are interpolated multilinearly. This accelerates the construction of nonlinear deformations and makes the method applicable for interactive applications. The method is compared to alternative approaches and the benefits are demonstrated on different kinds of input data.


Introduction
This paper is concerned with the construction of low-dimensional models describing the variability of a non-rigid shape from a set of example poses. Such models are used for templates-based surface reconstruction from incomplete data, data-driven shape editing, character animation, deformation transfer, and shape reconstruction from 2-dimensional data.
One approach for constructing such a model is to apply Principal Component Analysis (PCA) to the vertex positions of the example poses. The resulting PCA modes span a linear subspace of deformations that approximate the examples. Still, this is not an adequate shape model since the linear structure of the PCA cannot accurately represent shape deformations, which are usually highly nonlinear in the vertex positions. Furthermore, shapes do not deform under rigid transformations of R 3 , which can also not be captured by PCA. We consider shape manifolds that provide a shape representation invariant to rigid transformations. To this end, a point in a shape manifold corresponds to an equivalence class of poses, which differ only in rigid deformations. A Principal Geodesic Analysis (PGA) on such shape manifolds [HZRS18] can be used to obtain nonlinear models of shape variability that are invariant to rigid transformations of the example poses. A property of PCA and PGA is that they result in modes that involve all variables. In particular, they are not expected to be localized or decoupled. An alternative is the Sparse Localized Deformation Components (SPLOCS), introduced in [NVW*13]. The approach modifies the optimization problem underlying PCA by adding a sparsity inducing term to the objective. The resulting modes encode localized deformations that are intuitively meaningful.
In this paper, we introduce a Sparse Principal Geodesic Analysis (SPGA), which aims at combining the advantages of PGA with those of the sparse modes (cf . Figure 1). We consider the Riemannian manifold of discrete shells, which is a shape manifold of meshes equipped with a physically-based metric. The meshes are treated as the midsurfaces of thin shells and the metric measures the membrane and bending distortions of the shells. Our approach constructs a submanifold of the shell manifold and is invariant to rigid transformations of the example poses. Analogous to SPLOCS, we use a discrete L 1 term to enforce sparsity of the modes. However, in the nonlinear setting, the sparsity term has a different effect than in the linear setting. Since the metric measures membrane and bending distortions of the shells, the sparsity term forces the modes to describe localized deformations. For example, sparse modes can describe the movement of a human shoulder. Although the mode moves the whole arm, the distortions are active around the shoulder, hence, the distortion is localized (cf . Figure 2). In contrast, SPLOCS would not consider such a mode as sparse, because the sparsity refers to the local displacement instead of distortions. This property makes SPGA applicable to articulated shapes, which is a limitation of SPLOCS stated in [NVW*13]. Our approach builds on Nonlinear Rotation-Invariant Coordinates (NRIC), which represent a mesh by the vector that stacks the edge lengths and dihedral angles of all edges. The manifold of discrete shells corresponds to the set of such vectors, from which a mesh can be reconstructed. The SPGA modes are tangent vectors to the shell manifold at the mean shape. By applying the Riemannian exponential map to the space spanned by these modes, one obtains the submanifold of the shell manifold which approximates the input data.
The SPGA modes represent infinitesimal characteristic articulations of the mean and usually can be decomposed into lowdimensional almost completely decoupled subspaces, which is not possible for PGA modes. For example, for human models, one expects distinct, localized modes for the bending of elbow or knee whereas some more modes are required to represent shoulder ar-ticulation. By grouping the modes via their coupling and a simple summation of nonlinear coordinates corresponding to these groups, a practical approximation of the submanifold can be defined. Then, in an offline phase, one precomputes the exponential map on lattices in the low-dimensional decoupled subspaces. For given weights on the principal modes, they are then used in an online phase together with a multilinear interpolation and a very small number of Gauß-Newton iterations for the back projection onto the manifold. This yields a highly efficient discrete shell reconstruction and enables near real-time shape animation.
Contribution. We introduce an SPGA in the Riemannian manifold of discrete shells that combines the benefits of the PGA [HZRS18] and SPLOCS [NVW*13] in one model. Our main contributions are: • We use NRIC for PGA in the shell manifold resulting in a compact model for shape variability that is invariant to rigid body motions of the example poses and can handle nonlinear shape deformation. • By using a sparsity inducing term which reflects the shell metric, modes that describe localized deformations are obtained. This makes the analysis applicable to articulated shapes. • To solve the SPGA problem, we use a scheme that iteratively solves quadratic problems. In our experiments, this scheme outperforms the ADMM schemes used in prior work. • Based on the decoupling property of the sparse modes, we introduce a fast approximation scheme for the construction of a nonlinear deformation via superposition of deformations resulting from the decoupled groups of modes. This makes it possible to use our approach for interactive applications, which is not feasible for PGA.
The SPGA is evaluated on different data sets and compared to alternative approaches.
Organization. In Section 2, we discuss related work. The background on Nonlinear Rotation-Invariant Coordinates, geodesic calculus, and its variational discretization is revisited in Section 3. Based on this, in Section 4 the discrete PGA is introduced. Then, Section 5 presents the SPGA and Section 6 its computation via an efficient quadratic matrix programming approach. The superposition of large deformation submanifolds generated via grouping the SPGA modes is discussed in Section 7, whereas the lookup table based fast deformation synthesis is presented in Section 8. Experimental results are given in Section 9 and conclusions are finally drawn in Section 10.

Related work
Rotation invariant coordinates. For various tasks in geometry processing, it is beneficial to switch from the usual mesh representation in terms of vertex positions to an alternative representation that is adapted to the structure of the problem considered. An alternative to working with frame differences is to use the lengths of all edges and the dihedral angles between adjacent triangles [WDAH10]. In analogy to the classical differential geometry in which a surface is described by the first fundamental form, representing the metric of the surface, and the second fundamental form, characterizing curvature; edge lengths describe the metric of the mesh and the angles describe its curvature. Moreover, two meshes have identical edge lengths and angles if and only if they differ by a rigid transformation. Therefore, we call this representation the Nonlinear Rotation-Invariant Coordinates (NRIC). Integrability conditions, which guarantee that there are vertex positions realizing a given list of edge lengths and angles have been established in [WLT12]. The set of integrable NRICs forms a manifold and the tangential spaces at points in the NRIC manifold can be computed from the derivatives of the integrability conditions [SHHR20].
Since vertex positions do not exist for all combinations of edge lengths and angles, the question of how to project given lengths and angles to integrable NRIC coordinates arises. In [FB11] discrete shells, a model of elastic thin shells [GHDS03], is used for this purpose. Considering the mesh as the midsurface of a thin shell, membrane and bending distortions can be defined based on the edge lengths and angles of the mesh. To determine a mesh for a given list of edge lengths and angles, these are used to specify a rest configuration of the shell. The mesh is determined as a minimizer of the resulting discrete shell energy in the absence of constraints and external loads.
Shape spaces. Riemannian shape spaces are shape manifolds equipped with a Riemannian metric. In such spaces, concepts from Riemannian geometry can be used to derive tools for shape analysis and processing. This idea has been applied in computer vision, computational anatomy, and medical imaging and geometry processing. For a general overview of the topic, we refer to the textbook of Younes [You10]. A Riemannian metric on the space of triangle meshes and algorithms for computing the shortest geodesics, the Riemannian exponential, and parallel transport in this space were introduced and used for deformation transfer, shape interpolation and extrapolation in [KMP07]. A physically-based Riemannian metric was defined in [HRWW12; HRS*14]. The meshes are regarded to be discrete shells and the metric measures the viscous dissipation required to deform the shells.
Statistics in shape space. Principal Geodesic Analysis (PGA) [FLPJ04] is a generalization of PCA to data given in a Riemannian manifold. First, the Riemannian center of mass [Kar77] of the data points is computed and the data points are mapped to the tangent space at the center using the logarithmic map. Then a PCA in the tangent space is computed. Finally, the exponential map is used to map superpositions of the PCA modes back to the manifold. To perform nonlinear statistics for the analysis of shape variability, a PGA on the Riemannian manifold of discrete shells was introduced in [HZRS18]. In this paper, we extend this work by formulating the Shell PGA in the NRIC coordinates. These coordinates are well-suited for this purpose because the invariance under rigid transformations matches the structure of the Shell PGA.
Sparse PCA. PCA allows for analyzing and compressing sets of meshes representing animation sequences [AM00]. Varimax rotations can be used to alter the basis of the space spanned by the lowest principal components. For animation sequences, Varimax rotations were applied to obtain sparse components that are meaningful and could be used to support animators [MA07]. In [NVW*13], a sparsity inducing term is included in the computation of the principal components resulting in Sparse Localized Deformation Components (SPLOCS). Since the SPLOCS span linear spaces, larger deformations are not supported. To counteract linearization artifacts, a warping strategy, which represents the components using the deformation gradient and then treats rotational components nonlinearly, was introduced in [HYZ*14]. SPLOCS can also be formulated in the NRIC representation of meshes [WLZH17; LLW*19], which results in an improved analysis for articulated shapes. In contrast to our approach, the NRIC representation is treated as a linear space. We demonstrate that the PGA, which uses the NRIC manifold, its tangent spaces, the physically-based metric, and the Riemannian exponential and logarithmic maps, leads to improved results.
Data-driven shape editing. Statistics of shape variation are useful for a variety of applications in geometry processing and related fields. In this paper, we focus on data-driven shape editing.
Elasticity-based energies, such as As-Rigid-As-Possible [SA07] or discrete shells [GHDS03], are commonly used for mesh editing [SA07; HSvTP11; JBK*12]. These measure the energy stored in a deformation of a rest configuration. The choice of a rest configuration has a major influence on the deformation behavior, still, this choice is ambiguous. Data-driven shape editing uses a variable rest configuration, which is built from example poses, instead of a fixed one. In [FB11] the NRIC coordinates of example poses are linearly interpolated to obtain a model of a variable rest configuration for the discrete shells energy. While this works for a low number of examples poses, statistical shape models, such as a PCA in appropriate coordinates [GLL*16; WLZH17; LLW*19], PGA [HZRS18], or SPLOCS [GLY*20], are needed when dealing with larger numbers of example poses. Our experiments demonstrate the advantages of the proposed model for data-driven shape editing.

Background
We start with the necessary preliminaries to introduce our SPGA. First, we review NRIC as a useful representation for the geometry of a triangular surface to define the shape manifold of interest. Afterwards, we discuss the Riemannian structure and geodesic calculus on this manifold.
Nonlinear Rotation-Invariant Coordinates. In this paragraph, we briefly review the work by [WLT12] who introduced a discrete version of the fundamental theorem of surfaces. We consider a simply connected, triangular surface with a set of vertices V, edges E ⊂ V × V, and faces F ⊂ V × V × V. For nodal positions X ∈ R 3|V| , we denote by l(X) = (le(X)) e∈E the vector of edge lengths and by θ(X) = (θe(X)) e∈E the vector of dihedral angles.
To ensure that z = (l, θ) ∈ R 2|E| corresponds to a twodimensional triangular surface immersed in R 3 , it has to fulfill two admissibility conditions. The first condition is the triangle inequality, i.e. T f (l) > 0 for all f ∈ F where T f (l) = l i + l j − l k l i − l j + l k −l i + l j + l k for a face f ∈ F with edge lengths l i , l j , l k and the above inequality is to be understood componentwise. Furthermore, traversing from face to face along the fan of triangles at each vertex v in the set of interior vertices V 0 , one can reconstruct from l and θ the geometry of the fan up to rigid body motions. From this, one obtains a closing condition for the fan, which can be expressed as the admissibility condition Qv(l, θ) = 0, where Qv can be robustly and efficiently computed using quaternions. For details on this, we refer to [SHHR20].
Then, the manifold of all z ∈ R 2|E| which correspond to immersed triangular surfaces is given by where we collect all constraints in vector-valued functionals T = (T f ) f ∈F and Q = (Qv) v∈V0 . Following [SHHR20], we call M the NRIC manifold (Nonlinear Rotation-Invariant Coordinates). For z ∈ M, its tangent space is given by where DQ(z) ∈ R 3|V0|×2|E| is the Jacobian of Q. The conditions can be extended to higher-genus surfaces by including integrability conditions along non-contractible paths that generate the fundamental group, see also [WLT12].
The NRIC z enable a local description of shell deformations based on the variation of edge lengths encoding membrane distortions and the variation of dihedral angles encoding bending distortions. Frequently, in applications, these distortions and thus also the associated modulation of z components are localized. For example, on human models, in the vicinity of joints and active muscles which is in striking difference to vertex-based mesh representations (cf . Figure 2). Geodesic calculus. As in [HRWW12] we equip the manifold M with a Riemannian metric g reflecting the physical dissipation caused by the infinitesimal variation of a triangular surface interpreted as discrete shell. Following Rayleigh's paradigm, the metric acting on strain rates is related to the Hessian of the elastic energy, i.e. gz : is the elastic energy of a deformation of the discrete shell with coordinates z into the discrete shell with coordinatesz. It is composed of a membrane and a bending contribution with Wbend[z,z] = ∑ e∈E (θe −θe) 2 d −1 e l 2 e , where de = 1 3 (a f + a f ) for the two faces f and f adjacent to e ∈ E (a f is the area of f ) and Furthermore, µ and λ are positive material constants and G[z,z] is the Cauchy-Green strain tensor of the deformation, which is a function of the edge lengths on each face. The logarithmic term in the energy density Wmem ensures that the triangle inequalities remain fulfilled for finite-energy deformations.
For the Riemannian distance, we have that dist 2 (z(0), z(1)) = E and straightforwardly ). Thus, we replace the squared Riemannian distance by the computationally much cheaper elastic deformation energy and obtain the discrete path energy for discrete paths (z 0 , . . . , z K ) ∈ (R 2|E| ) K+1 . Then, discrete Kgeodesics are defined as minimizers of the discrete path energy for given end points z 0 and z K and the associated system of Euler-Lagrange equations is given by

Principal Geodesic Analysis
In this section, we will discuss the generalization of Principal Component Analysis on linear spaces to Riemannian manifolds. In this context, the Riemannian logarithm and exponential map enable to transform nonlinear, large variations on the manifold to infinitesimal variations given as tangent vectors in a particular tangent space and vice versa. Given the already introduced discrete geodesic calculus, we use discrete counterparts of the continuous logarithm and the continuous exponential map. We define the discrete logarithm at some discrete shell z 0 as is a discrete initial velocity of the discrete geodesic (z 0 , z 1 , . . . , z K ) based on time step size τ = 1 K . Here, PT z0 M is the orthogonal projection onto the tangent space Tz 0 M with respect to the metric gz 0 , which is extended to R 2|E| per its definition. The corresponding exponential map Exp z0 (v 0 ) for v 0 ∈ Tz 0 M, as the inverse mapping, computes a discrete geodesic (z 0 , z 1 , . . . , z K ) for given discrete, initial velocity v 0 . In fact, we first compute z 1 for given z 0 and v 0 via orthogonal projection of z 0 + v0 K onto M with respect to gz 0 , i.e. z 1 := The projection PM can be computed using a Gauß-Newton scheme [FB11]. Then we recursively compute z k+1 solving After K − 1 steps this yields the requested discrete geodesic and the discrete exponential, i.e.
For given points z 1 , . . . , z N ∈ M, the continuous Fréchet meanz is defined as the arg min z∈M ∑ n=1,...,N dist 2 (z n , z). In what follows, we replace continuous geodesics used in the definition of dist(·, ·) by discrete K-geodesics to define a discrete (Fréchet) K-meanz. Of course, these definitions depend on the choice of K, however, we always used K = 8 for Exp and Log and did not observe any benefit by increasing it in our experiments, see also Section 9. Concerning the mean, we obtained already very good results for the 1-meanz, which can be considered as the elastic mean [RW09].
Now, we are in the position to introduce the discrete Principal Geodesic Analysis (PGA) [FLPJ04] for M. Again, let z 1 , . . . , z N ∈ M be a given data base of discrete shells.
Discrete PGA algorithm 1. compute the meanz of the input shapes z 1 , . . . , z N . 2. compute v n = Logz(z n ) for n = 1, . . . , N. 3. perform a PCA on the tangent space TzM with gz as the underlying scalar product to compute the dominant J principal modes of variation u 1 , . . . , u J ∈ TzM.
The actual PCA computation (3.) can be formulated as where V ∈ R 2|E|×N is the matrix containing the v i as columns and the columns of U are the modes u 1 , . . . , u J . The approximation error is measured in a Frobenius norm · g weighted according to gz. Taking into account the minimizing property and the fact that v n ∈ TzM, one easily checks that all u j are in TzM.
For each u j , the curve t → Expz(t u j ) can be interpreted as the jth nonlinear mode of shape variation of the input data. Moreover, considering the whole subspace U = span{u 1 , . . . , u J } of the tangent space TzM the PGA allows us to define a submanifold of the NRIC manifold M which approximates the input data.
In comparison to linear PCA in R 2|E| , i.e. using the (projected) linear mean and linear differences, the nonlinearity allows us to better capture the articulation of the input shapes. This means for the same linear approximation quality, i.e.
, the actual approximation of the input shapes z n is more accurate using nonlinear PGA and the exponential map than using linear PCA and orthog- onal projection on M. A comparison of this nonlinear PGA and a PCA in NRIC coordinates is given in Figure 3.
In general, Expz(t u j +s u i ) = Expz(t u j )+Expz(s u i ) by the nonlinearity of the exponential map. Nevertheless, after replacing standard PGA modes by sparse modes in NRIC coordinates in the next Section 5, the resulting manifold M[U] itself will be approximated by linear superposition of very low-dimensional submanifolds in Section 7 -thus allowing for an efficient deformation synthesis in Section 8.

Sparse Principal Geodesic Analysis
By definition, conventional PGA outputs an orthonormal set of dominant deformation modes. These modes usually have large support as already the orthogonality implies largely overlapping supports of pairs of modes to ensure the annihilation in the gz product. Motivated by the quest for more efficient parametrizations of data approximating submanifolds mentioned above and by the intuition that, in many applications, we are interested in sparse, spatially localized dominant modes with widely disjoint supports. For example, on human shapes, it seems natural to deal with dominant deformation modes which are supported only on a foot, the surrounding on a knee, or the area of a shoulder. The choice of NRIC, which represent local distortions via local coordinate chances, enables such a sparsity on articulated shapes in our approach.
To define sparse deformation modes, we follow the common approach of considering a sparsity-inducing regularization term. However, using an unweighted l 1 -norm on tangent vectors would be problematic as they consist of edge lengths and dihedral angles varying on vastly different scales. To this end, we introduce a family of weighted L p -norms Working without the tangent space constraint leads to modes with a mismatch between the length and the angle component. The dihedral angle support of a sparse linear PCA mode in R 2|E| coordinates which has no length component at all is shown on the left. This leads to mesh artifacts already in case of short time extrapolation as seen in the middle picture. Furthermore, this is prohibitive also for the approximation of shapes, which can be seen in the last picture which shows the failure to approximate an input shape, which was easily handled even by the linear PCA in R 2|E| , even though the linear approximation quality is about 96%.
wherede is the area associated with e in the meanz. The scalar product associated with the L 2 -norm has already been used in [FB11] to define a quadratic energy on lengths and angles. To simplify the notation in the following, we collect the weights in a vector m ∈ R 2|E| such that u 1 = |m u| 1 , where means entry-wise multiplication. With this L 1 -norm at hand, we generalize the PGA to a Sparse PGA (SPGA) as follows.
We dropped the orthonormality constraints as we do not aim for an orthonormal basis anymore but favor sparsity. To ensure that the magnitude of the coordinates of the u j does not simply shrink to achieve a small L 1 -norm while the weights grow accordingly, we introduce a bound on the magnitude of the weights' entries. Different from standard PGA, we have to impose that the modes u j are indeed tangent vectors in TzM, which ensures that they yield admissible, infinitesimal deformations of the meanz. Dropping this constraint DQ(z)u j = 0, one would leave the NRIC manifold M and there is no guarantee that u j represents a geometrically admissible variation of the underlying triangular mesh. To illustrate this, we added in Figure 4 just the L 1 -regularization to the linear PCA. As we can see, this led to components of edge length and dihedral angle variation with non-matching support. Their usability for approximation and the synthesis of new shapes is quite limited. In fact, it prohibits the use of extrapolation via the exponential map.
Finally, we can again by virtue of the exponential map define a submanifold of the NRIC manifold M containing nonlinear deformations.

Quadratic Matrix Programming
We solve (4) by alternatingly solving for U while keeping W fixed and vice versa. When U is fixed, solving for W becomes a straightforward quadratic optimization problem. However, solving for U is more difficult due to the regularization term and the tangent space constraint.
To solve for the sparse modes efficiently, we rephrase (4) for fixed weights W as a quadratic matrix programming problem [Bec07]. To this end, we first reformulate the nonlinear matrix optimization problem, and then relax the L 1 -penalty term as in [Tib96; BH17] to obtain a quadratic problem.
Matrix optimization. Let us begin by reformulating the matrix optimization problem to prepare the relaxation. We denote by G ∈ R 2|E|×2|E| the matrix representation of the metric gz , i.e. G = W ,22 [z,z]. By the definition of the (weighted) Frobenius norm, To rephrase the regularization term, we denote by M ∈ R J×2|E| the matrix whose rows are all equal to the weight vector m. This allows the following reformulation We express the tangent space constraint using the quaternion integrability operator Q, i.e.
Together, this yields for given W the problem While this problem already closely resembles a QMP problem, the choice of the L 1 -term still renders it nonlinear.
Solving the matrix optimization problem by relaxation. Next, we relax the L 1 -penalty (cf. [BH17;Tib96]) to obtain the quadratic problem with inequality constraints. To this end, we introduce non- This way, using the componentwise triangle inequality |U| ≤ U + + U − we obtain the estimate for the L 1 -term. Analogously, adapting the other terms in (5) we finally arrive at the relaxed quadratic matrix programming problem The above triangle inequality is an equality, i.e. |U| = U + +U − if U + U − = 0. Hence, the solutions U + ,U − of (6) have disjoint support because, otherwise, we could move values from one matrix to the other without changing the value of the objective term while decreasing the regularization term. Thus, U = U + − U − is indeed the solution of the original problem. To solve the relaxed problem, one can first vectorize it (cf . [Bec07]) and then utilize off-the-shelf software for quadratic programming.

Effective Submanifold Approximation
In many applications, the sparse modes (e.g. see Figure 5) can at least partially be split into subsets with pairwise decoupled support. Recall, for human models one expects distinct, localized articulations for certain joints and muscle groups.  Figure 6 and Figure 10).
Hence, the supports of corresponding groups of sparse modes representing these articulations are well separated. Furthermore, as pointed out in Figure 7, the support of these tangent vectors in TzM ⊂ R 2|E| is only very moderately extended under application of the exponential Expz even for large K. Given such a strong separation of supports of two subspaces Ua and U b one observes that for ua ∈ Ua and u b ∈ U b . This is certainly not the case for subspaces corresponding to modes with overlapping supports, which is illustrated in Figure 6 using the inset modes. In fact, the defect of a potentially non exact separation of supports after applying the exponential map can be measured using where (·, ·) 2 is the scalar product associated with the weighted L 2norm introduced in Section 5, see also Figure 7. Before, we have already argued that the support of the tangent vectors is only moderately extended under the exponential map. The same then holds for the coupling of pairs of tangent vectors, which allows us to measure the separation of modes and induced submanifolds Expz(Ua) and Expz(U b )using which can also be seen when comparing Figure 8 and Figure 7.
This observation motivates the following splitting approach. Let us suppose that the space U ⊂ TzM spanned by the sparse principal modes u 1 , . . . , u J can be written as the direct sum of L subspaces U l for l = 1, . . . , L, i.e.
where each of these subspaces is spanned by a subset of the principal modes. In fact, as discussed above, we suppose that due to our sparse PGA approach the subspaces U j can be chosen to have at least approximately pairwise disjoint supports. For each vector v ∈ U, we have a corresponding decomposition v = v 1 + v 2 + . . . v L with v l ∈ U l . Now, we consider a superposition of exponential maps via summation of NRIC coordinates

Efficient Deformation Synthesis
So far, we introduced an approximation of the submanifold by splitting the application of the exponential map into applications on smaller subspaces. It remains to show how these exponential maps and the needed projection can be approximated efficiently to enable the fast synthesis of nonlinearly deformed shapes on the submanifold. We propose an approach based on sampling in the preprocessing phase and multilinear interpolation in the online phase.
To this end, we first consider a single subspace U l of dimension d with basis u l 1 , . . . , u l d . In this subspace, we take a lattice {∑ d n=1 αn τ u l n | α ∈ Z d } with mesh size τ ∈ R. For α ∈ Z d , we denote by v α = ∑ d n=1 αn τ u l n the corresponding lattice point. Now, an arbitrary v = ∑ d n=1 wn(v) u l n ∈ U l lies in a cell of the lattice with ∈ Z d elementwise. In two dimensions, this would mean that α(v) identifies the lattice point to the lower-left of v. To approximate Expz(v), we consider a piecewise multilinear interpolation Z τ,l of the values at lattice points, i.e.
where the w β (v) are the multilinear interpolation weights (cf . to cheaply evaluate the piecewise multilinear interpolation in the online phase.
Considering the whole subspace U = U 1 ⊕ . . . ⊕ U L , we use for general v ∈ U the splitting v = v 1 + . . . + v L and repeat the interpolation above for each subspace. Then, we evaluate the superposition of these interpolations The necessity of the projection is demonstrated in Figure 9.
Surely, when dealing with the submanifolds Expz(U l ) in our computational setup, we are interested only in compact subsets of the subspace U l representing plausible deformations. Furthermore, in the applications, it turned out to be sufficient to restrict to subspaces U l of low dimensionality between one and four. Hence, the multilinear interpolation of precomputed samples on a rectangular grid becomes a feasible option. A comparison of the exact submanifold M[U l ] = Expz(U l ) and the approximation M τ [U l ] for a single two dimensional subspace U l is shown in Figure 10. The efficiency of this approach based on multilinear interpolation of precomputed samples is evaluated in the next section and in particular in Table 4.  Figure 10: Approximation of the exponential map by multilinear interpolation of lattice points for two strongly coupled modes (cf. Figure 5). On the left, we show the lattice (grey) in the twodimensional subspace along with randomly sampled points colored according to the logarithm of their relative approximation error in the energy, i.e. W[Expz(v), zτ[v]]/W[z, Expz(v)] . On the right, selected extrapolated lattice points are shown in purple and we highlight the quality of the approximation of shapes (purple point clouds) by our deformation synthesis approach (yellow) for samples with relatively high approximation error. The linear approximation quality is measured as V − UW 2 / V 2 and weighted L 1 -norm and overlap are computed as described before. Noticeably larger values for λ led to vanishing modes. In this case, we chose J = 10 and λ = 250.

Experimental Results
We have implemented our method in C++ with the Geometric Optimization And Simulation Toolbox (GOAST) [HS*20], where we use the Eigen library [GJ*10] for numerical linear algebra and CHOLMOD [CDHR08] and UMFPACK [Dav04] from the SuiteSparse collection as direct linear solvers. For the quadratic problem (6), we use MOSEK [MOS20] as efficient off-the-shelf solver.
We use a multi-resolution approach to enable fast computations while retaining the possibility to produce high-quality deformations. In this approach, the input shapes are coarsened by one of two methods: For the human model and the horse input data, we used an iterative edge collapse approach based on minimizing the quadric error metric [GH97] computed in groups [MG03] to preserve the dense correspondence between input shapes. For the hand and face input data, we remeshed a reference input shape using OpenFlipper [MK12] and computed coarse representations of the other input shapes with the same approach as in the prolongation described next. In both cases, the coarse results were prolongated to the fine level using a representation of the fine mesh vertices in terms of intrinsic positions and normal displacement with respect to the coarse mesh and computed on a reference shape similar to [KMP07]. See Table 4 for a comparison of original and coarse resolutions of the different input data.
Selection of number of modes and sparsity. In the computation of modes, we have two parameters which are currently chosen manually based on heuristics: The number of modes J and the sparsity weight λ in (4). This is linked to a trade-off between the approximation quality of the model, the sparsity and overlap of support of the modes, and the size of the model. In Figure 11, we compare these quantitatively for a number of choices and also to PGA. For the approximation, not only the average overlap is important but also the resulting subspace dimensions. Thus, a higher number of modes is not necessarily beneficial.
Selection of coupling spaces. The selection of the subspace decomposition U = U 1 ⊕ · · · ⊕ U L , i.e. the grouping of the sparse tangential modes based on their coupling, is an important step to obtaining high-quality results with our proposed method as exemplarily pointed out in Figure 6. In our implementation, we employ spectral clustering [SM00] of a set of modes using the sparse tangential coupling D(u i , u j ) as the underlying similarity measure. To this end, we use the coupling matrix of the sparse modes, cf . Figure 8, as affinity or similarity matrix. The number of clusters is then determined such that the resulting dimensionality of the subspaces U l is in the range 1, . . . , 4, see also Figure 12, and thus small enough for the lattice generation described in Section 8 to be computationally feasible. The decomposition into subspaces of different dimensions used in our examples is listed in Table 1. The impact of the coupling of sparse tangential modes on the decoupling of the associated nonlinear submanifolds is depicted in Figure 7. Examples of groups of modes extracted from the SCAPE data set are shown in Figure 1. The examples illustrate that modes that move the same part of the body, for example the left leg in the leftmost groups, can be in different groups. The reason is that the distortions induced by the modes are located in different areas of the body, in the hip and knee regions for the leftmost groups.  for increasing γ, each again solved with Gauß-Newton. We choose the starting value for γ and increase it by a fixed factor until the mesh fulfills the handle positions up to a pre-defined tolerance.
In Figure 14, we compare this mesh editing approach to other state-of-the-art, data-driven approaches. To compare our approach to [LLW*19], we replicated their approach without localization, i.e. using a spatially constant l 1 -term in the mode computation, and with a fixed penalty parameter (see the right-most example in the top row). The other two methods, [FB11] and [GLL*16], are applied as described in the original publications. Our method delivers plausible deformations on-par with [HZRS18], while being nearly two orders of magnitude faster. Some details are even superior, for example, the bending of the right arm, shown in the close-ups, appears more natural in our results with the sharp elbow contour being better preserved.
Using the superposition of nonlinear subspaces represented by zτ instead of the linear combinations of modes in R 2|E| allows our method to generalize better for edits requiring deformations far from the input data set. For example, in Figure 15, we consider the challenging example of dragging the right hind leg outwards of a horse shape, while the input data only consists of a galloping sequence. Here, besides the penalty approach used in Figure 14 we used hard constraints for the vertex positions at the handles. We compare the editing via a linear basis and hard constraints in R 2|E| (left), which leads to unnatural bending similar to the one reported PGA SPGA Example Off-Online Off-Online Fig. 13 (Interpolation) 3min 3s 90min 60ms Fig. 14 (Fitting) 2min 919s * 660min 10s (100ms) Table 2: Comparison of runtimes for PGA and SPGA. In the interpolation case, we report the average time to evaluate the spline in Figure 13 at 120 evenly spaced points. In the fitting case, the total amount of time to produce the results in Figure 14 is reported and for SPGA, in brackets, also the time for incremental editing steps. * as reported in [HZRS18] 2 1 2 2 2 3 2 4 2 5 2 6 10 −2 10 −1 Relative error Log 2 1 2 2 2 3 2 4 2 5 2 6 K 10 −2 10 −1 Exp 2 1 2 2 2 3 2 4 2 5 2 6 2 7 10 −2 10 −1 Exp•Log−Id Figure 16: Convergence of time-discrete exponential and logarithm for K → ∞ on the SCAPE dataset. For the logarithm and exponential, we show the relative error in the weighted L 2 -norm using Kmax = 128 as pseudo ground truth. Furthermore, due to the tangent space projection the discrete exponential does not coincide with the inverse of the discrete logarithm. But, we observe that it is an approximation with decreasing relative error in the weighted L 2 -norm for K → ∞.
in [LLW*19, Figure 17], with editing based on the large deformation interpolation zτ[·], both for hard constraints (middle), and using the quadratic penalty method (7) (right), where the latter two yield far more natural results.
Comparison PGA to SPGA. In comparison to PGA [HZRS18], our SPGA computes sparse modes at the expense of approximation accurateness. Hence, to achieve the same approximation quality with SPGA as with PGA one will typically need more modes (cf . Figure 11). The sparsity, however, is crucial for the submanifold approximation in Section 7, as otherwise the resulting error would be large (cf . Figure 6). When working with PGA, we cannot use our efficient deformation synthesis and instead have to directly use the nonlinear exponential map. In editing applications, when additional handle constraints come into play, this implies that we have to evaluate derivatives of the discrete exponential and hence use methods from PDE-constrained optimization, see [HZRS18]. This results in online runtimes which are far from interactive rates ( Table 2). In contrast, our approximation built on SPGA allows for interactive rates but requires more runtime in the offline phase for the quadratic optimization and sampling of the exponential map.
Comparison to ADMM. Another method commonly used to compute sparse deformation components is the alternating direction method of multipliers (ADMM), cf . [NVW*13; HYZ*14; WLZH17]. We have also investigated this approach to solve for the sparse modes. To deal with the tangent space constraint in (4) one needs to add a third term to the commonly used ADMM problem representing this constraint as a convex indicator function.   Table 4: Overview of used datasets with corresponding preprocessing and online timings listing from left to right: the data set, the number of input shapes N, the number of vertices of the full mesh |V| and the coarse mesh |Vcoarse|, the number of modes J, the timings in the preprocessing phase to compute the meanz, a discrete logarithm Logz, one solution of the quadratic matrix programming problem, a discrete exponential map Expz, and the timings in the online phase to compute the multilinear interpolations Zτ, a Gauß-Newton iteration used for the projection on the manifold M, and a Gauß-Newton iteration required in the mesh editing context. * obtained from six shapes by computing weighted elastic averages as usual for ADMM one alternatingly evaluates the proximal mappings for each term. A more detailed explanation of this can, for example, be found in [BPC*11]. However, we observed a slower convergence, i.e. compared to our approach outlined in Section 5 the objective value was larger after the same number of outer iterations while still requiring more time per outer iteration. We provide objective values and runtimes for the two methods in Table 3.
Time discretization. To use the Riemannian exponential map and logarithm as tools in our method, we needed to discretize them in time, cf . Section 3. The discretization is controlled by the number of steps K in a discrete geodesic and converges for K → ∞ to the corresponding continuous notions, which is shown in [RW15]. We demonstrate this convergence empirically for the SCAPE dataset in Figure 16 and also saw that for K = 8 we already achieve an acceptable error. Thus, we chose it for all experiments reported above.
Timings. In Table 4, we list detailed timings of all components of our method. For the quadratic problem, we report runtimes with enabled parallelization, while it was disabled for all other timings. The offline phase was performed on a workstation with two AMD EPYC 7601 CPUs, the high number of cores allows to efficiently parallelize the offline phase by computing many logarithms or exponential maps in parallel. All other results were computed on a laptop with an Intel Core i7-9750H CPU.
In the preprocessing phase, we used 25 outer iterations, i.e. solving once for the weights and modes, in each of the examples. Precomputing the exponential for lattice points (cf . Section 8) required about 3000 (horse), 4000 (hands), 20000 (human), and 30000 (face) evaluations for the respective examples. In the online phase, evaluating our deformation synthesis zτ for a random v ∈ U using the mean shape as initialization needed in all examples at most ten Gauß-Newton iterations. For larger, instant changes, such as in Figures 14 and 15, we need multiple increase of the penalty parameter in (7) and thus a larger number of iterations. Concretely, the computation for the result in Figure 14 took about ten seconds while for Figure 15 it took about five seconds. In applications with incremental changes, we only need a very small number of Gauß-Newton iterations. For example, each step of a spline for the shapes in Figure 13 evaluated at 120 evenly spaced points required two to three iterations, while for the incremental handle editing (shown in Figures 1 and 17) even one was sufficient. This makes the method applicable for interactive applications.

Conclusions
We introduced a framework for data-driven surface processing that allows synthesizing large, nonlinear deformations based on a Sparse Principal Geodesic Analysis on the space of discrete shells and enables close to real-time articulation of surface models and interactive mesh editing. Thereby, we use the Nonlinear Rotation-Invariant Coordinates to deal with the space of discrete shells as a Riemannian submanifold of R 2|E| . The approach consists of an offline and an online phase. -In the offline phase, we first use a discrete geodesic calculus to compute a mean shape and the discrete Riemannian logarithms of the input surfaces. From these logarithms, we identify sparse, dominant modes of shape variability. Second, the modes are grouped with respect to their overlap and by superposition of the discrete Riemannian exponential map acting on the span of each group a submanifold of the shape manifold, which fits the input data, is defined. Third, we compute a lookup table of sample shapes extrapolated this way. -In the online phase, the synthesis of large nonlinear deformations is performed by summing of multilinear interpolations of samples from the lookup tables and finally projecting these close approximations by very few steps of a Gauß-Newton scheme onto the shape manifold. We tested our framework on different shape data sets, evaluated the performance quantitatively and compared it to different state of the art methods. Limitations and challenges. The current implementation can only handle simply connected surfaces and input data which is in dense correspondence. A generalization to higher-genus surfaces would require the handling of additional, nonlocal integrability conditions. Computing the geodesics requires sufficient mesh quality of the input shapes to prevent local minima. We witnessed on the SCAPE dataset, that our method sometimes got stuck in such local minima. This, however, did not degrade the performance of the method due to its rareness. Still, extending the geodesic calculus to less quality meshes would increase the method's usability to a wider variety of datasets. The method still requires a projection onto the shape manifold after the superposition of groupwise extrapolations, constructing an approximation omitting this would improve performance and enable more flexible applications. Furthermore, a multiscale approach tailored to NRIC could help to further increase the method's performance. We have reported error measurements for the superposition of extrapolated groups of modes. Here, an a priori control of the error would enable a theoretically rigorous grouping of modes. Moreover, it would be interesting to combine the variational sparsity approach based on L 1 -norms with a variational mode untangling as in [SBBG11]. Finally, recent hierarchical tensor product approximation tools could allow generalizing the currently single scale grouping of modes and support an improved and more efficient manifold approximation.