{'fraction_non_alphanumeric': 0.03017454596848279, 'fraction_numerical': 0.007278637552116459, 'mean_word_length': 5.0283280085197015, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Regularized Harmonic Surface Deformation

Harmonic surface deformation is a well-known geometric modeling method that creates plausible deformations in an interactive manner. However, this method is susceptible to artifacts, in particular close to the deformation handles. These artifacts often correlate with strong gradients of the deformation energy.In this work, we propose a novel formulation of harmonic surface deformation, which incorporates a regularization of the deformation energy. To do so, we build on and extend a recently introduced generic linear regularization approach. It can be expressed as a change of norm for the linear optimization problem, i.e., the regularization is baked into the optimization. This minimizes the implementation complexity and has only a small impact on runtime. Our results show that a moderate use of regularization suppresses many deformation artifacts common to the well-known harmonic surface deformation method, without introducing new artifacts.


Introduction
Surface deformation is an important task in geometry processing. Deforming models involves interactive modeling sessions driven by a user, who deforms an object by manipulating a subset of the surface vertices. Linear deformation methods [BS08] have proven effective in this context as they often create plausible and realistic-looking deformations, while still allowing for interactive runtimes. Deformations are usually modeled as the minimizers of specific deformation energies that are defined locally at each point of the domain and measure a specific distortion property.
Harmonic surface deformation is a well-known linear gradient-domain method introduced by Zayer et al. [ZRKS05] that uses a differential surface representation to perform global mesh deformations. Deformation constraints are smoothly distributed over the entire mesh using harmonic functions and surface details are preserved in the reconstructed deformations. Furthermore, the method is parameter-free.
Linear methods such as harmonic deformations are prone to artifacts, as linear energy terms are inad-equate to accurately model the non-linear physical forces and processes involved in a deformation. Common deformation artifacts include flipped triangles (in the case of planar surfaces), protruding triangles, degenerate elements, volume loss as well as local and global shape distortion. Many of those artifacts occur close to the deformation handles. Figure 1 shows an example.
A number of non-linear correction methods [Lip12,AL13,SKPSH13,KABL14] exists that suppress these artifacts in planar or volumetric settings. These methods are very powerful as they guarantee artifactfree deformations, but they typically have a big impact on the runtime. Most importantly, they are not applicable to deformations of surfaces. Martinez Esturo et al. [MRT14] introduce an alternative linear energy regularization method: a quadratic regularization term is proposed that is strongly coupled to the problem-specific deformation energy. For a number of problems, this regularization yields to artifact-free results, albeit it cannot be guaranteed. Technically, this energy regularization requires only minor modifications to the algorithm with little impact on the runtime. The amount of regularization can be adjusted using a single parameter.
In this work, we apply linear energy regularization to the harmonic surface deformation method of Zayer et al. [ZRKS05]. Hereby, we follow the general ideas of Martinez Esturo et al. [MRT14]. We demonstrate that energy regularization enhances harmonic deformation results and suppresses a variety of artifacts. Our main contributions are: • We provide an energy-regularized formulation of harmonic surface deformation.
• We refine the discretization of the energy differential operator of Martinez Esturo et al. [MRT14] for better estimates in high curvature regions.
• We evaluate the effectiveness of our approach. In particular, we demonstrate that moderate use of energy regularization improves deformation results by resolving artifacts without introducing new ones.
This paper is structured as follows: we discuss related shape editing techniques and correction methods (Section 2) and review harmonic deformations as well as energy regularization (Section 3).  The base of the model is fixed, and each sphere represents a deformation handle, which was rotated. Harmonic surface deformation is susceptible to protruding triangles and surface self-intersection near small deformation handles. Low amounts of regularization (β) suppress these artifacts.
introduce our approach to linear energy regularization for harmonic surface deformation (Section 4). We perform a qualitative analysis of our results (Section 5), followed by quantitative evaluation and discussion (Section 6). Lastly, we present our conclusions and outlook for future work (Section 7).

Related Work
The goal of interactive surface deformation is to create meaningful deformations while preserving surface properties such as local details and curvature. Linear deformation methods play a major role in this area, since they provide the interactivity and often produce plausible deformations. Most often, linear methods represent the surface using its differential properties [Sor06]. One can distinguish these methods with respect to their sensitivity regarding rotation and translation. Linear techniques often cannot guarantee that the used deformations are smooth everywhere [JBPS11]. This leads to deformation artifacts such as flipped triangles, protruding elements and volume loss due to rotations.
Artifacts can be avoided by improving the smoothness of the transformation interpolation field. However, bi-or tri-harmonic weights create additional local extrema in the interpolation field that lead to unintuitive deformations results. Jacobson et al. [JWS12] tackle this problem by forcing a desired topology for the interpolation field. This requires solving a nonlinear conic problem.
Serveral correction methods have been explored as another means to reducing artifacts in various geometry processing tasks. Lipman [Lip12] presents a generic tool for constructing orientation preserving (i.e., no triangle flips allowed) triangle mesh mappings, while limiting worst-case conformal distortion. This method has non-interactive run times and is only defined for planar meshes. Schüller et al. [SKPSH13] propose a specialized optimization based on a barrier energy function to repress zero-area elements and flipped triangles at interactive rates. Their iterative scheme solves for an injective mapping to a new mesh configuration. They guarantee inversion-free mappings of planar triangular and volumetric tetrahedral meshes. Aigerman and Lipman [AL13] extend [Lip12] to volumetric meshes. Their algorithm takes a deformation created by common deformation techniques and returns a similar deformation that is injective and minimizes the distortion of the mesh volumetric elements. The method is not interactive. Most recently, Kovalsky et al. [KABL14] present a method based on linear matrix inequalities for restricting the range of singular values. It enables, e.g., bounded distortion mappings of planar or volumetric domains, but is also computationally to expensive for interactive applications.
These correction methods guarantee that deformations are inversion-free, and in some cases even protrusion-free. However, they all require solving non-linear systems, which results in loss of interac- tivity for moderate to large meshes. In contrast, we follow the recent linear approach to regularization by Martinez Esturo et al. [MRT14], which has no significant impact on runtime. While we cannot guarantee artifacts-free deformations, our method successfully suppresses usual deformation artifacts. Furthermore, the correction methods mentioned above are not applicable to surface meshes embedded in R 3 . In contrast, our method is well-define for surface meshes.

Background
In this Section, we continue to review the formal details required in our work. We consider triangulated surface meshes M = (T , V, E) defined by sets of vertices i ∈ V, oriented edges E ⊂ V 2 , and triangles T ⊂ V 3 . Coordinates of vertices i ∈ V are denoted by x i ∈ R 3 . A missing subscript either indicates a vector of stacked coefficients, e.g., x ∈ R 3|V| , the vector of stacked vertex coordinates x i , or a matrix of component-wise coefficients, e.g., X ∈ R |V|×3 . For a triangle t ∈ T , X t ∈ R 3×3 denotes the column-wise concatenation of the coefficients of its vertices. Using the notations we denote (squared) vector and matrix norms that are induced by symmetric and positive definite matrices N. (Tr(·) denotes the trace of a matrix.) For the piecewise linear functions on M a discrete gradient operator G ∈ R 3|T |×|V| can be assembled from local per-triangle gradient operators G t : for triangles t = (i, j, k) ∈ T with normalized normals n t , the local gradient operators are given by see, e.g., [BS08]. Given a scalar function on M defined by the vertex-based coefficients u ∈ R |V| , G u is the vector of stacked and constant per-triangle gradients. Note that gradients computed by G are defined in a common coordinate system.

Harmonic Guidance for Surface Deformation
Zayer et al. [ZRKS05] propose a variant of gradientdomain deformations in which local deformation constraints are propagated using harmonic functions. Deformed surfaces are reconstructed from manipulated surface gradients by minimizing the global deformation energy subject so suitable boundary constraints. Here, A t denotes the area of triangle t, ||M|| 2 F = Tr M T M is the (squared) Frobenius norm of M, and Z t ∈ R 3×3 are prescribed component-wise guidance gradients that are constant per triangle. Dense guidance gradients are computed from user-specified transformations associated to a set of handle regions. In [ZRKS05], harmonic functions h(x) given by solutions of the Poisson equation L h = 0 are used for the global propagation of the sparse set of given transformations. Here, L = G T A G is a discretization of the Laplace-Beltrami operator [BS08], in which A is a diagonal matrix of replicated triangle areas. Please, see Zayer et al. [ZRKS05] for further details on the quaternionbased propagation of transformations using harmonic functions. As minimizers of (3) are characterized by the same differential operator, in practice a single factorization of L can be used to perform both transformation propagation and energy minimization.
Drawbacks. Linear deformation methods are susceptible to various artifacts. This is because the physical energies involved in deformations are nonlinear by nature, and are only approximated by linear methods [BS08]. Specifically, harmonic surface deformation is susceptible various deformation artifacts, which are in particular related to the size of the handle regions. Small deformation handles are likely to cause protruding and intruding triangles. Large deformation handles cause local shape distortion on the boundary between the constrained and free mesh areas. Other artifacts include volume loss close to deformation handles when the deformation includes large rotations and local surface intersection. It was observed [MRT14] that these artifacts correlate with large spatial variation in the optimized energy on the domain.

Linear Energy Regularization
Martinez Esturo et al. [MRT14] propose a generic linear energy regularization scheme that suppresses geometric artifacts in a number of different applications. It is applicable to regularize problem-specific squared energies of the general form Here, E ∈ R n|T |×d|V| is a problem-specific linear energy operator that maps unknown functions u to triangle-constant local energies of dimension n, c are problem-specific energy constants, and A is a diagonal matrix of replicated triangle areas that performs domain-wide integration of triangle-constant quantities. Energies E P are regularized by introducing a regularization term that measures squared variations of local energies.
For piecewise constant local energies, pointwise energy variations are estimated by the sparse differential operator D. For each pair of neighboring triangles, it can be discretized along all internal edges e ∈ E i from the set of nonboundary edges E i ⊆ E: let l(e) and r(e) denote the left and right triangle at e, respectively. Then, for scalar local energies (n = 1), the nonzero coefficients of D are given by , for all internal edges e ∈ E i and triangles t ∈ T . For vector-valued local energies (n > 1), the differential operator is given by a component-wise replication, which can be expressed as D⊗I n using the Kronecker product ⊗ and the n×n identity matrix I n . The constant estimates of pointwise local energy variations are integrated using the diagonal matrix B of replicated internal edge lengths.
The total regularized energy is given by a weighted combination of both terms that can be expressed compactly using the β-weighted norm The amount of regularization is steered by β ∈ [0, 1).
Note that this formulation of energy regularization is also valid for energies in the components of the unknown functions u, which are then given by E β (U) = ||E U − C|| 2 W β . Please, see Martinez Esturo et al. [MRT14] for further details and applications on this energy regularization scheme.

Enhancing Harmonic Surface Regularization
We continue to show that the concept of energy regularization is applicable to harmonic surface deformation in a straightforward way. For this, we rewrite the component-wise deformation energy (3) to the equivalent formulation using the global gradient operator G, the matrix Z of all stacked prescribed gradients, and diagonal matrix A of replicated triangle areas. Note that the local energies correspond to the to the summed terms of (3). Comparing our problem-specific energy (8) and the regularizable generic energy (4), we obtain the correspondences that the generic energy operator E is given by the gradient operator G, the constant energy term C is given by the gradient field Z, and the dimension of the local energies is n = 3. Hence, the energy-regularized version of the deformation energy (8) is given by by simply applying the weighted norm W β for energy integration and smoothness estimation. Technically, to apply regularization, this substitution of A for W β allows for a straightforward implementation. Specifically, the remainder of the original harmonic surface deformation approach is unchanged, in particular the harmonic function-based transformation propagation for the guidance gradients Z.
Our experiments demonstrate that this simple energy modification suppresses a variety of deformation artifacts of the original energy formulation (see Section 5). Still, the original discretization of the energy differential operator D is defined independently of the local surface curvature, which leads to poor estimates of energy variation in high curvature regions. We continue to provide a refined differential operator discretization that is based on surface curvature and yields better estimates of energy variation.
Curvature-based Energy Differential Operator. Martinez Esturo et al. [MRT14] estimate the variation of local energies by finite differences of their respective local energy residuals (Section 3.2). In our application of gradient-domain deformations, the energy residual E t of a triangle t is given by the deviation between the deformed mesh gradients and the gradients of the guidance field: The (squared) energy variation between two neighboring triangles t 1 and t 2 is now simply given by ||E t1 − E t2 || 2 F . This estimation works well if both triangles are aligned and therefore also the residuals in the columns of E t1 and E t2 live in the same tangent space, e.g., for planar meshes. However, if the tangent spaces are not aligned, this estimation is likely to break. The effect is particularly noticeable at sharp edges of a curved surface, where areas of low energy values along the edge are surrounded by higher energy values. This results in a less smooth energy distribution across edges of high curvature.
We adjust the discretization of the operator D to also be applicable to high curvature regions by locally compensating for curvature. Similar local alignments of tangent spaces is used, e.g., by Crane et al. [CDS10], for the computation of connections relating neighboring tangent spaces.

Re
For the n = 3 dimensional local energy residuals (10), our refined energy differential operator D R is a |E i | × |T | matrix of 3×3 block matrices. For each internal edge e ∈ E i , we use the rotation R e that aligns the normals of its left and right triangles n r , n l such that n r = R e n l . Then, the nonzero blocks of D R are given by for all internal edges e ∈ E i and triangles t ∈ T . This way, D R compensates for local curvatures: residuals that are similar relative to their respective tangent spaces are not estimated to be different anymore. Note that this operator refinement simplifies to the original formulation for planar meshes. It reduces the sparseness of the resulting linear system only by a constant factor, which does not impair resulting runtimes. However, it is only applicable to this particular case of n = 3 dimensional local energy residuals. For different values of n the discretization should therefore fall back to the D operator of [MRT14].
Implementation. The operators G, D R , A, and B as well as L are assembled once when the surface mesh is loaded. For each deformation of the model the guidance field Z is computed and the corresponding normal equations of (9) are solved for the coordinates X of the deformed mesh. The norm W β is assembled for a given β value using the refined differential operator D R . After elimination of positional hard boundary constraints, the system (11) is symmetric positive definite and it is solved using a Cholesky factorization with fill-in reducing reordering [GJ + 10]. Similar to [ZRKS05], this factorization can be reused for different guidance fields Z as long as the handle configuration and β values are unchanged.

Results
We evaluate our approach qualitatively on a number of different models. The models presented have 3 − 36.6k vertices. Deformations are performed using single and multiple handles, varying sizes of the handle regions and the regularization weight β.
Energy Visualization. We use color to visualize the squared magnitude of local triangle-constant energies ||E t || 2 F . Energies are linearly mapped to the color space (shown in Figure 1) by setting the maximum of the color interval to correspond to the 95th percentile of the energy values. The top 5% of the energy values is clipped to allow visualization of the variation of lower energies with higher contrast. resulting in protruding triangles and local surface intersections near the deformation handles. Regularization (β > 0) suppresses these artifacts. Similar artifacts to protruding triangles, albeit of smaller magnitude, occur if handle regions include up to several dozen vertices. In the Cactus 2 example shown in Figure 3, single vertices on the base and the top are fixed, and a single vertex on the side of the cactus is used as a deformation handle. The deformation suffers from large protruding triangles. Introducing regularization corrects this artifact. For both deformations, even low amounts of regularization result in a smoother energy distribution, affecting more triangles in the mesh. This way the optimization favors more global changes to the mesh instead of local concentrations of energies, which result in the observed local artifacts.
Large Deformation Handles. Harmonic surface deformation is commonly used with large deformation handles, as careful design of these constraints can lead to pleasing deformations. However, the boundary of large handle regions is also susceptible to artifacts. The Cactus 1 model in Figure 5 is based on a benchmark deformation from [BS08]. The base of the model is fixed, and the top is rotated and translated. Using the original method by Zayer et al. [ZRKS05], the boundary between the constrained and free mesh regions exhibits strong changes of the directions of mesh normals, distorting the local shape. Our regularization reduces distortions of local geometry close to this boundary, creating a smooth transition between the constrained and free mesh regions.
Strong Rotations. Strong rotations, especially on elongated limbs, are a challenge for a number of deformation methods [KS12]. Using harmonic surface deformation, strong rotations usually cause loss of volume. This effect is illustrated in Figures 4 and 2. The foot of the Horse is rotated. In the resulting deformation, most of the lower leg suffers from volume loss, creating a "candy wrapper"-like artifact. Low values of regularization are effective at preserving some of the limb's volume to create more plausible results. In the Cow deformation, the horns are rotated upwards while the front of the face is fixed. This deformation is even more challenging, as the region deformed is small and the mesh geometry is rather coarse. Regularization helps to restore some of the lost volume.

Evaluation and Discussion
In this section, we perfom a quantitative evaluation of our approach. Deformations are created for β ∈ [0, 1), using a higher sampling rate for lower β values, for which we usually observe the strongest changes of the deformation. Please, see the accompanying video for the deformations at these β values.

Maximal Isometric and Conformal Errors.
Two error measures proved to be most useful for evaluating deformation quality (see, e.g., [LZX + 08]): the local per triangle isometric error e iso t is given by the sum of squared deviations from 1 of the singular values of the deformation gradient. The local per triangle conformal error e conf t is computed by half of the squared sum of the pairwise differences between the singular values of the deformation gradient. The maximal isometric (conformal) error, max(e iso t ) (max(e conf t )), is given by the maximal value of the isometric (conformal) error over all triangles. Both of these error measures indicate strong distortions of the mesh geometry. They are also loosely related to artifacts such as protruding or intruding elements, degenerate triangles, and surface self-intersections. We also examined the total (integrated and normalized) isometric and conformal errors as possible indicators for deformation quality evaluation, but these proved ineffective. Figure 6 shows the behavior of the maximal isometric and conformal errors for different values of β. In almost all tested deformations, regularization strongly decreases the maximal isometric and conformal errors. The behavior of these error measurements is different for the Cactus 1 example of Figure 5, as they slightly increase with regularization. The reason is that deformations defined using large handle regions usually don't result in artifacts associated with large local isometric and conformal errors. Additionally, we note that no new local maximal errors are introduced for moderate β values, meaning that deformation quality does not derogate for higher β values. We also confirm this behavior in all other tested examples, as the total deformation energy of the regularized deformation defined in (3) stays within the same order of magnitude as the original deformation.
Space of Regularized Deformations. We observe that the initial introduction of regularization (β > 0) has strong effects on deformations suffering from strong artifacts for β = 0. After this initial reaction interval, regularization creates gradual changes to the mesh geometry, indicated both by our deformation results (e.g., in Figures 2 and 5) and by the behavior of the error measures max(e iso t ), max(e conf t ) in Figure 6. As the regularization energy becomes more dominant with increasing β values, max(e iso t ), max(e conf t ) also slightly increase. For very high values of β > 0.9, the regularized energy formulation strongly deviates from the original problem, which can create new artifacts. In our experience, choosing a fixed value of β ≈ 0.2 usually suppresses artifacts effectively without negatively affecting the mesh geometry. This means that a constant regularization can simply be added to existing implementations without exposing users to a new parameter. In addition, although the linear system (11) has to be refactored when β changes, examining different β values can usually be done at interactive rates due to the high performance of sparse linear solvers. Hence, the space of regularized deformations can also be explored interactively.
Curvature-based Differential Operator. Figure 7 demonstrates the benefits of using our curvature-based energy differential operator on the Accordion mesh with highly curved edges. The method by Zayer et al. [ZRKS05] suffers from local shape distortion near the deformation handle. The differential operator by Martinez Esturo et al. [MRT14] estimates the energy variation between neighboring triangles ineffectively, resulting in global shape distortions. Our curvature-based operator estimates the energy differential more effectively, resulting in a deformation which has the same geometry as the solution by Zayer et al. [ZRKS05], but suppresses the local artifact. The deviation of Martinez Esturo et al. [MRT14] from our results increases with higher regularization weight. For smooth meshes, using our curvature-based differential operator has negligible effects, as vertex coordinates were only affected marginally.
Performance. Regularization affects the runtime in two ways: it reduces the sparsity of the linear system (11) being solved and it requires a one-time operation to setup of the curvature-based differential operator D R . Our measurements on an Intel Core i7 2.2GHz system indicates that the effect of regularization on interactive runtime is insignificant for all tested meshes: For example, factorization time of the linear system for the Cactus model (|V| = 10k) is ≈ 0.04 seconds for both β = 0 and β > 0, its solution time is ≈ 0.003 seconds. Similarly, factorization time for the Hand model (|V| = 36.6k) is ≈ 1.13 seconds for β = 0 and ≈ 1.14 seconds for β > 0 with a solution time of ≈ 0.1 seconds. Hence, changing the regularization weight to examine different regularization weight can be done interactively.

Conclusions
In this work, we have provided an energy regularized formulation of harmonic surface deformation. Our approach expands the capabilities of the original method, allowing the creation of artifact-free deformations for a wider range of deformation constraints and handle configurations. Our formulation of a curvature-based differential energy operator improves the estimation of energy differentials in highcurvature mesh regions. This reduces geometric dis- Harmonic surface regularization [ZRKS05] creates an artifact near the deformation handles. Martinez Esturo et al. [MRT14] suppress the artifact, but distorts the global mesh shape. Our curvature-based differential operator creates a smooth, artifact-free deformation without distorting the mesh geometry.
tortions introduced by the original energy differential operator around these regions. The evaluation of our results demonstrates that even low regularization weights can effectively suppress many deformation artifacts without negatively affecting the performance of the original method. In addition, no new artifacts are created, Future work. An interesting direction for future work is the application of energy regularization to other (non-linear) surface deformation approaches, e.g., to the work of Jacobson et al.
[JBK + 12], who observe similar artifacts. In addition, since artifacts are usually localized close to handle regions, a local energy smoothness formulation could achieve better control on deformation artifacts. Hand Horse Cactus 1 Cactus 2 Figure 6: Graphs of the maximal isometric and conformal errors over the regularization weight β. Regularization typically decreases these errors drastically already for small values of β. Note the logarithmic scale on the y-axes.