ARAP Revisited Discretizing the Elastic Energy using Intrinsic Voronoi Cells

As‐rigid‐as‐possible (ARAP) surface modelling is widely used for interactive deformation of triangle meshes. We show that ARAP can be interpreted as minimizing a discretization of an elastic energy based on non‐conforming elements defined over dual orthogonal cells of the mesh. Using the intrinsic Voronoi cells rather than an orthogonal dual of the extrinsic mesh guarantees that the energy is non‐negative over each cell. We represent the intrinsic Delaunay edges extrinsically as polylines over the mesh, encoded in barycentric coordinates relative to the mesh vertices. This modification of the original ARAP energy, which we term iARAP, remedies problems stemming from non‐Delaunay edges in the original approach. Unlike the spokes‐and‐rims version of the ARAP approach it is less susceptible to the triangulation of the surface. We provide examples of deformations generated with iARAP and contrast them with other versions of ARAP. We also discuss the properties of the Laplace‐Beltrami operator implicitly introduced with the new discretization.


Introduction
As-rigid-as-possible (ARAP) surface modelling [SA07] is a popular and practically relevant approach for interactive deformation of triangle meshes with significant impact on research in geometry processing [LZX*08, BKP*10, TCL*13, BDS*12].The basic approach is to minimize a deformation energy subject to a set of constrained vertices.The deformation energy is governed by penalizing the deviation of a vertex star from transforming rigidly: This energy is minimized over the set of vertex positions V and the per-vertex rigid transformations R i simultaneously.Deviation from rigidity is measured by comparing the rotated original edge vectors êij to the edge vectors e i j in the deformed mesh among the edges incident on a vertex.The weights w i j control the influence of each edge.Sorkine and Alexa suggest to use the cotan weights known from discrete Laplace operators for triangle meshes [PP93,MDSB03].Their argument is heuristic and based on the observation that the choice of diagonals in a rectangular grid should not affect the deformation (cf. Figure 1).The fact that the weights are negative for non-Delaunay edges is a well known problem, because it may cause the per-vertex energy term to become negative.
This may explain why it has become more common to use a modified version of the ARAP energy, following a proper discretization of the continuous deformation energy suggested by Chao et al. [CPSS10].The idea of this energy, similar to ARAP, is to penalize the deviation of the gradient (coordinate wise, i.e.Jacobian) of a deformation mapping f : R 3 → R 3 from a rigid transformation.The smooth version of this energy can be evaluated exactly in closed form for a linear deformation function and constant rigid transformation over a triangle [PP93]: (i, j)∈t cot α t i j e i j − Rê i j 2 . (2) Here, α t i j is the interior angle in triangle t opposite edge (i, j), as in the definition of the cotan Laplace operator.
Modifying the ARAP energy to sum over the triangles incident on vertex i leads to the so-called spokes-and-rims version (Figure 2)  [SA07,LSHXY22].The constrained vertices (coloured in red) are displaced orthogonal to the plane.The results exhibit the asymmetry in the energy of the spokesand-rims approach.As the initial mesh is Delaunay the iARAP approach coincides with the original ARAP approach.where in addition to the edges incident on vertex i also the edges in the link (i.e. the boundary of the vertex star) contribute.While the cotan weights still may be negative for non-Delaunay edges, the energy exactly expresses the non-negative contribution of each triangle and cannot become negative in any part.
The addition of the rim edges introduces a potentially detrimental asymmetry in the cells that are governed by a common rigid transformation (see Figure 1).In Section 3, we interpret the different versions of ARAP as different discretizations of the elastic deformation energy by Chao et al. [CPSS10].This shows that the spokesand-rims version of ARAP is based on the mesh triangles, while the original ARAP may be considered using orthogonal dual cells.If the mesh is Delaunay the orthogonal dual cells are intrinsic Voronoi regions and this works well.Non-Delaunay edges, however, lead to dual cells that are not properly immersed.We ensure valid dual cells by considering the intrinsic Delaunay triangulation, whose dual is the intrinsic Voronoi diagram and the dual cells are properly immersed.This construction leads to the iARAP energy.The necessary computation of local rotations requires that we need to make the intrinsic triangulation available extrinsically.In Section 4, we explain how to represent the intrinsic Delaunay edges as polylines living on the original triangle mesh using barycentric coordinates.
The modified ARAP discretization differs from the original version in that it is guaranteed to be non-negative for each vertex star.This remedies the artefacts observed in the original ARAP version.It behaves similar to the spokes-and-rims variant, yet avoids unwanted asymmetries for meshes whose discretization fails to reflect the symmetries.We provide visual results and quantitative data in Section 5.
A particular appeal of ARAP surface modelling is a straightforward block-descent scheme for minimizing the energy: Rotations are computed locally by non-linear optimization, vertex positions are computed globally by solving a linear system with fixed system matrix.The system matrix is a discrete Laplace operator, for both the original and spokes-and-rims version of ARAP.The matrix resulting from our approach is different.In other words, the original ARAP and the spokes-and-rims variant have the same left-hand side in the linear system, and the problems resulting from the possibly negative energy are resolved by modifying the right-hand side of the system by spokes-and-rims.In the iARAP approach, the right-hand side is similar in spirit to the original ARAP, but the left-hand side is modified.We analyse the properties of the left hand side as a discrete Laplace operator in the spirit of the discussion by Wardetzky et al. [WMKG07] in Section 6.As it turns out, the new construction allows trading symmetry for a maximum principle.
Lastly, we briefly discuss how our approach relates to other modifications of the energy, in particular those that address the dependence of bending on the discretization.We also provide an outlook on using the implicitly introduced Laplacian operator for other applications.

Background and Notation
We consider the immersion of a fixed triangle mesh M that is being deformed by a function f : M → R 3 .The immersion M is defined by the positions vi of the vertices V and the assumption that triangle i jk ∈ T is realized as the convex hull of its vertices.
The deformed mesh M is the result of applying f to the vertex positions, that is v i = f(v i ), and then realizing the triangles as the convex hull of their mapped vertices, identical to the definition of M. Note that we explicitly avoid asking that f also maps the interiors of the triangles to the convex hull of the mapped vertices.This gives us the freedom to consider functions f that are not necessarily piecewise linear or continuous on the given mesh.
Given a realization of the mesh, edge vectors can be computed as e i j = v j − v i .This can be done for both, the original e i j and the deformed version of the edge êij .The circumcentre of triangle i jk in the original mesh is ĉijk .The dual cell i of vertex i is a polygon formed by the circumcentres of the triangles incident on i. See Figure 3 for an illustration.The dual cell is immersed if all primal edges i j are Delaunay; otherwise, it is self-intersecting, with dual edges êij corresponding to non-Delaunay edges having opposite orientation.
Discrete laplace operators.It will be useful to recall the constructions of the cotan Laplace operator.There are two substantially different constructions leading to the same result.They have recently been contrasted by Alexa el al. [AHKSH20] in the context of pointing out that the constructions differ for tetrahedralizations.The primal discretization is based on linear basis functions associated with the triangles of the mesh.We provide slightly more detail on the dual approach for context of our derivation of the ARAP energy; for more detail we ask the reader to consult the extensive literature on the subject and also see the derivation in Section 6.
Assuming the mesh consists of Delaunay edges, consider the dual cell i.The integrated Laplacian over i turns into a sum over triangles spanned by vertex i and consecutive circumcentres, e.g.v i , c i jk , c i jk .For each triangle we need to make an assumption of the function.The common approach is to consider only the (values in the) vertices along the primal edge, for example i and j.Other vertices, such as k and k , only contribute to the triangles associate with corresponding primal edges.This means, the function value assumed in c i jk generally differs depending on the primal edge being considered.Therefore the underlying finite element space is nonconforming: The function is discontinuous along the edges of the triangles comprising the dual cell.In contrast, the finite element space for the primal construction is based on linear elements on triangles that agree on edges and is conforming.In general there is no reason for two different discretizations to agree and the fact that that the primal and dual constructions do lead to the same discrete operator for the Laplacian on triangle meshes is a coincidence.More generally, there are no inherent advantages or disadvantages of conforming or non-conforming finite element spaces and both are commonly used in different domains.

Intrinsic triangulations.
The dual edge e i j of non-Delaunay edges points in the "opposite" direction and leads to negative edge "lengths".This results in negative weights for the cotan Lapace-Beltrami operator and may even result in negative cell areas.This can cause problems when the cotan discretization is used [SA07,CWW17].The intrinsic Delaunay triangulation of the mesh yields a principled solution to this problem.An intrinsic edge between two vertices on a polyhedron is a geodesic edge on the piecewise linear surface.In the realization of the surface it may be considered a polyline (compare Figure 7).Bobenko and Springborn [BS07] define the Laplace-Beltrami operator based on intrinsic Delaunay edges.They show that any polyhedron has a unique intrinsic Delaunay triangulation that can be constructed by intrinsically flipping non-Delaunay edges.Other constructions of the intrinsic Delaunay triangulation are possible [LFXH17] and they differ in their worstcase time complexity.In practice the computation is observed to be linear in the number of elements of the mesh for real-world models [SSC19,LFXH17].For representation there exist several data structures [FSBS06, SSC19, GSC21] with support for different demands in applications.

Derivation of the ARAP Energies
We start from the elastic deformation energy introduced by Chao et al. [CPSS10]: The main idea of ARAP is to restrict the rotations R. It seems that all published variants of ARAP restrict R to be piecewise constant.The variants differ in the regions over which R is considered constant, and in the restrictions imposed on the deformation function f. 3) is

As shown by Chao et al. [CPSS10], a necessary condition for minimizers of the energy in Equation (
which can be thought of as analogous to the situation for minimizers of Dirichlet energy having vanishing Laplacian.This is, in fact, the global step common to different versions of ARAP, and may be used as an explanation why all of them are based on the cotan Laplace-Beltrami operator on the LHS in the linear system.
The simplest version of ARAP results from using constant rotations per triangle.This choice has the drawback that bending across edges is not penalized, but it may be useful when other constraints are considered simultaneously [LG15,ZG18,LJ21].In this case, R is computed in the local step based on the three vertices of a triangle, and divergence is computed in the usual way [PP03,Hir03], considering the columns of the rotation matrices as vector-valued functions for each of the three coordinates.
Explaining the original ARAP energy and the spokes-and-rims version in this framework requires making somewhat "unusual" choices.Both are based on fixing a rotation for each vertex and it seems natural that the boundary of this region intersects each edge at its midpoint which we define for both the original and the deformed geometry.In this way we have set m i j to be the midpoint of the mapped vertices v i , v j , meaning we assume that f maps edge midpoints to edge midpoints.
With this assumption, we can interpret the rotations of spokesand-rims ARAP as being defined over the polygon formed by the edge midpoints.While it was originally suggested to compute R i from the vertex star of vertex i, using the polygon of edge midpoints is identical, except for an irrelevant scale factor.This means, we can interpret spokes-and-rims ARAP as considering the energy over only 3/4 of each triangle as depicted in Figure 4. Looking at derivations for spokes-and-rims ARAP, we see that divergence of  the rotations in vertex i is computed by first averaging the three rotations over each triangle, and then computing the divergence from the triangle-wise quantities as usual.
Using the polygon of edge midpoints makes the discretization dependent on not only the positions of vertices, but also the connectivity of the mesh.A natural alternative, avoiding at least the dependence of the choice of edges, would be to use the Voronoi regions of the vertices.We depict the area considered for the two cases in Figure 5. Indeed, it is possible to interpret the original version of ARAP in this way.Note that the Voronoi cell around vertex i is spanned by a set of circumcentres, if the mesh is Delaunay.Our assumption will be that the deformation f may map these circumcentres freely and also differently depending on what triangle incident on the circumcentre we consider.In other words, we allow f to be discontinuous.Consider the primal edge (i, j) and assume it has the Delaunay property.The part of the dual (Voronoi) cell i corresponding to the edge is the triangle vi , ĉijk , ĉijk .This triangle is intrinsically flat, but in its realization it is comprised of two triangles, connected along the edge vi , mij .Let us focus on the triangle vi , ĉijk , mij .By our choice of f, it will be mapped to the triangle v i , c i jk , m i j .Note that c i jk is the result of mapping ĉijk considering the sub-triangle incident on vertex i and edge (i, j).When considering the sub-triangle on the same edge but incident on vertex j the edge is considered in the opposite direction and the mapped circumcentre is c jik .Likewise, sub-triangles incident on edges ( j, k) and (k, i) lead to four more circumcentres c jki , c k ji and c ki j , c ik j .Because we allow f to be discontinuous, all six mappings may be chosen different.This degree of freedom can be exploited in order to minimize the deformation energy.This means, applying f to the piecewise linear surface M results in decomposing each triangle into six triangles, not necessarily connected across edges, illustrated in Figure 6.Note that although f is discontinuous our definition of M means that its triangles are flat and different from the original triangles mapped by f (except for the edges).
Equation (2) describes the contribution of this smaller triangle to the deformation energy.Notice that edge c i jk − v i is opposite to a right angle, so its contribution is zero.The edge mij − ĉijk transforms to m i j − c i jk .For any choice of R i minimizing the energy, we can set c i jk so that R i ( mij − ĉijk ) = m i j − c i jk .In this way we have chosen f so that the contribution of this edge to the energy is minimal, namely zero.This leaves only the transformed edge v i − m i j contributing to the energy.This is true for both sub-triangles of vi , ĉ i jk , ĉijk .The overall contribution of this triangle to the energy then is: Here we have exploited the fact that the angles ∠v i c i jk m i j and ∠v i v k v j are identical (∠v i c i jk m i j is half the central angle ∠v i c i jk v j , which is twice the inscribed angle ∠v i v k v j ).Considering also the triangle v i , m i j , c i jk we get the suggested choice for the weights w i j in Equation 1.
Lastly, for computing divergence of R consistent with the original ARAP formulation we consider the diamonds (v i ĉijk v j ĉijk ) around vertex i (instead of the triangles).Each diamond consists of the surface spanned by two primal vertices incident on an edge and the corresponding two dual vertices (circumcentres).The rotation per diamond is computed by averaging the two rotations associated with the primal vertices.Then divergence is computed over the diamonds, similar to computing it over triangles [PP03].
Note that this derivation is based on the assumption that dual cells are immersed and that basing the realization on circumcentres yields Voronoi cells.If the mesh is not Delaunay, these assumptions are invalidated.Our idea to overcome this problem is to use the intrinsic Delaunay triangulation to define intrinsic Voronoi cells as the elements used for the discretization.This, in contrast to the approach of Chao et al. [CPSS10], keeps the spirit of computing divergence from the diamonds (albeit of the intrinsic Delaunay triangulation) for the RHS in Equation 4, but it changes the Laplace operator on the LHS.

Intrinsic ARAP Energy
The intrinsic Delaunay triangulation has been used mostly in a context where it is sufficient to compute the length of the edges [CWW17].In order to measure the deformation of the surface we need to consider the realization of the intrinsic Delaunay edges.Note that each intrinsic Delaunay edge (i, j) is a polyline living on the piecewise linear surface.We denote the line segments of this polyline as s t i j , where t denotes the triangle of the original mesh that carries the segment.In case the intrinsic edge coincides with an edge of the mesh the choice of t is not unique, but this has no consequences.The intersections p uv i j of the polyline with the extrinsic edge (u, v) are the start and end points of a segment: For an illustration see Figure 7.
Note that this notation is not well defined, as a single intrinsic edge may intersect the same extrinsic edge multiple times [SSC19].
To avoid clutter we omit this complication in our presentation -it has no impact on the calculations.
We want that the energy depends on the vertex positions in the deformed state (and not the segment vectors).Note that we can express the intersections of an intrinsic edge with a mesh edge in terms of barycentric coordinates relative to the mesh vertices.A segment is the difference of two such barycentric coordinates, living on the same triangle.This means we can write for the segments in the rest state where V = (v 0 v1 . ..) is the matrix that stores the vertex positions in the rest state column wise and b t i j is a sparse vector representing the segment in barycentric coordinates.The fact that b t i j encodes a vector (and not a point) means its elements sum to zero instead of one; and since both endpoints are on the same triangle it has at most three non-zero elements.If the segment coincides with a mesh edge it has only two non-zero elements.In particular, this is how original mesh edges are represented if they are Delaunay.
Similar to fixing the weights w i j based on the rest state geometry we also fix the barycentric representations b t i j .In this way we now simply replace the straight segment e i j in the ARAP energy with the corresponding set of line segments {s t i j }: The scalar factors σ t i j account for the fact that the intrinsic edge has been divided into smaller segments.We determine these factors by asking that the energy contribution of an edge is independent of its subdivision into segments (we drop subscripts referring to the edge to avoid clutter).Consider the intrinsic view, in which the edge is straight, we have s t = s t e e and note that t s t = e .Equality of energy for the edge and the sum of the segments yields as required.We further define w t i j = σ t i j w i j .

Derivatives
For the optimization it is necessary to calculate derivatives.The gradient of the energy with respect to the vertex position v k is where (b t i j ) k is the k-th element of b t i j , as s t i j = Vb t i j .
As (b t i j ) k is only non-zero for segments within the star of v k and each segment s t i j also appears negated in the energy, as segment s t ji , with the same weight w t i j = w t ji , we can write the derivative as where S k defines the set of directed segments within the extrinsic star of k, not the intrinsic.

Optimization
We can use the same block-descent optimization method as the original ARAP method [SA07].The optimization is divided into a local step that optimizes {R i } given V and a global linear step optimizing V given {R i }.
For given vertex positions, the local step consists of calculating optimal rotations R i for each cell centred around ṽi .Our approach differs from the original version only by the definition of the edges associated with cell of i.We arrive at the following covariance matrix: Optimal rotations R i can be computed in different ways -we use the polar decomposition The global step is directly given by setting the gradient w.r.t.vertex positions (Equation 14) to zero: Note that this is the same linear system of equations for every dimension and the system matrix only depends on the initial mesh.
Recall that we can write s t i j as Vb t i j .This suggests that the system matrix is where B contains all the barycentric coordinate vectors b t i j of all segments and D w is a diagonal matrix with the weights corresponding to the segments in B. This representation is similar to the definition of the polygon Laplacian by Alexa and Wardetzky [AW11, Equation 3] for Delaunay meshes, as B reduces to the co-boundary operator in this case.
The matrix L is symmetric by construction and positive semidefinite because the weights in D w are positive.So we can use the Cholesky factorization similar to the other ARAP methods.

Results
We implemented the iARAP method using Eigen [GJ*10] and libigl The iARAP discretization requires computing the intrinsic Delaunay triangulation in a preprocess, adding to the time required for precomputation (Table 1).After this, each step in the optimization is comparable to the original ARAP and spokes-and-rims ARAP.Only the calculations of the segments comprising an intrinsic edge are more involved.We store the precalculated sparse barycentre vectors b t i j and then the additional effort is a single sparse matrix-vector multiplication.The overhead grows with the number of intrinsic edge crossings as shown in Table 2.The resulting overhead is negligible in practice, see Table 1.
In Figure 8, we show the decay of the energy over the iterations.We observe that iARAP behaves similar to spokes-and-rims ARAP and slightly better than the original ARAP.Anecdotally, all methods seemed to offer the same degree of interactivity during modelling sessions.
The visual differences compared to the original ARAP on non-Delaunay meshes, as shown in Figure 9, are noticeable.Due to negative energy contributions the original ARAP exhibits creases, for example at the ears of the Bunny.In fat, the creases already appear  without displacing the constraints.As the spokes-and-rims ARAP and iARAP energies are well defined, we regularly observe smooth and similarly looking deformed meshes, see Figure 10.
The difference between spokes-and-rims ARAP and iARAP clearly shows on meshes consisting of nearly regular quads that have been triangulated by inserting diagonals.Figure 11 shows the deformation of a cube with the top of the cube moved backwards.From the side all three methods look the same, which is why we only show one instance.In the front view we see asymmetries in the spokesand-rims variant.This effect may also be observed on real world shapes, such as the balloon depicted in Figure 12.     [ZJ16].Notice the symmetric gap between the legs in the input.As the original ARAP energy diverges we stopped the optimization after a few steps (and the image shows only a subset of the resulting mesh).Spokes-and-rims ARAP and iARAP both converge, but only iARAP preserves the symmetry.

Figure 14: On the base of the lion, intrinsic edges (red) cross a lot of extrinsic edges (blue). While the results of spokes-and-rims ARAP and iARAP are roughly similar, spokes-and-rims slightly bends the base.
Asymmetries are also visible in symmetric shapes whose representation as a triangle mesh fails to represent the symmetries.The stool from thingi10K [ZJ16] has symmetric legs.Symmetric deformations should not remove this symmetry, yet only iARAP exhibits the desired behaviour (Figure 13).Subtle differences can be observed in quite generally when the triangulation of the surface is far from being Delaunay or, in other words, if the intrinsic Delaunay edges intersect many extrinsic edges.The base of the lion in Figure 14 is triangulated with many non-Delaunay edges.While iARAP is oblivious to this triangulation, spokes-and-rims ARAP slightly bends the flat region, presumably due to numerical issues.

Laplace-Beltrami Operator
Both the original ARAP method and spokes-and-rims ARAP yield the cotan Laplace-Beltrami operator as the left hand side of the lin- ear system of equations for the global optimization step.The iARAP discretization leads to a different matrix.We discuss this matrix interpreted as a Laplacian operator (Section 3) in terms of the list of properties introduced by Wardetzky et al. [WMKG07].
As mentioned before, Equation 17 shows that L is symmetric and positive semi-definite.
The operator is local in the sense that L i j = 0 if there is no edge between v i and v j .This is the case because the derivative with respect to v i only depends on the star of v i .On the other hand, the weights depend on the intrinsic Delaunay triangulation, which is global, so the operator is only weakly local in the sense that moving vertices may affect the weights of vertices that are far away [WMKG07,AHKSH20].
Linear precision is equivalent to (Lv) x = 0 for each interior vertex v x if all vertices are straight-line embedded into the plane.Given such a vertex v x we prove linear precision by first grouping all segments s t i j ∈ N x in the star of vertex x according to their intrinsic edge e i j .It is clear that each group defines a straight polyline starting and ending at the boundary of the star or starting/ending at v x .The two cases are visualized in Figure 15.We represent a polyline as an ordered list of extrinsic points (p 0 , p 1 , . ..) using sparse barycentric vectors (similar to segments, see Section 4): The contribution of a polyline P that starts and ends at the boundary of the star of v x is where p t s is the start point and p t e the end point of the segment s t i j .Since all points lie in a common plane, all segments of the polyline between i and j point in the same direction.Therefore is the same for all segments and equals This simplifies the contribution since successive segments share a point: where p t 0 s is the start and p tn e is the end vertex of the polyline.As both are on the boundary of the star of x this contribution is 0 as (b i ) x = 0 for all points i on the boundary of the 1-ring.
For the second case we only consider polylines P ∈ P x that start at v x , as the segments are directed and the opposite direction contributes the same way: (b t i j ) i w t i j s t i j .
Note that x = i and that all polylines in P i that start at i contain only one segment.This is clear as the end point of a segment in a triangle t that starts at v i must be on the edge opposite to v i .We can again rewrite the term due to the linearity of the vertex positions and the fact that the polylines start at v i as where N i describes the intrinsic Delaunay neighbourhood of i and p tn e is the end vertex of the polyline.As p tn e is on the boundary of the star of i by definition and (b i ) i = 1 we end up with: In the planar embedding this is the same contribution as the intrinsic Laplace-Beltrami operator and according to Stokes' theorem this adds up to zero [AHKSH20].Concluding, the iARAP operator yields the same result as the intrinsic cotan Laplace-Beltrami for all linear functions embedded in a plane, by weighting only the extrinsic neighbourhood.
From the results of Wardetzky et al. [WMKG07] it follows directly that the "positive" weight property cannot be satisfied or else we would have a contradiction with their main result.Note that the definition of the sign of the Laplacian varies in the literature.In our setup we would ask that all off-diagonal elements of L are negative.There is no reason why BB T should not contain positive and negative weights, as b t i j contains negative and positive values.An overview of the properties is shown in Table 3. Summarizing, the properties of the Laplacian implicitly defined by the intrinsic ARAP discretization of the elastic energy of Chao et al. [CPSS10] has properties that are identical to the cotan Laplacian.We want to stress, however, that the Laplacian corresponding to iARAP is generally not the cotan Laplacian.Experimentally we find that if both have (unwanted) negative off-diagonal entries, the iARAP discretization leads to smaller absolute values.This may help to alleviate unwanted effects of the wrong sign of these coefficients.
Alternative derivation.One may derive the operator implicitly defined by iARAP in a slightly more general way using Stokes's theorem: Here, f is an arbitrary function, the Laplace-Beltrami operator, ∇ the gradient operator, C a cell on the domain of f and n the outward pointing normal of the cell boundary δC.
In a discrete setting the function f is represented as vector f.If we define C to be the intrinsic Voronoi cell i, the dual edges e i j resemble the boundary.By assuming a constant gradient along the dual edge e i j , we can approximate the integrated gradient in normal direction: We arrive at the (integrated) intrinsic cotan Laplace-Beltrami operator: where N i is the intrinsic neighbourhood of i.
In the case that i has only knowledge of the function values of its extrinsic neighbours this approximation is not possible anymore.If we still want to resemble the intrinsic Voronoi cell, we could approximate the intrinsic edges by tracing outwards.The intersections of the extrinsic link with the intrinsic Delaunay edges could be used to approximate the function value of the intrinsic neighbour by linearly interpolating the two extrinsic neighbours on the link.
Using our notation for the barycentric coordinates of the traced edges we find the following term where N i is the intrinsic neighbourhood of i and b j p V is the intersection of the intrinsic edge (i, j) and the boundary of the vertex star of i.
In fact this is exactly the iARAP Laplace-Beltrami formulation considering only polylines that start at i.If we set the contribution to all other polylines within the star of i to 0 we arrive at an extrinsic approximation of the intrinsic Delaunay Laplace-Beltrami operator derived from finite volumes.This new Laplace-Beltrami operator L + has positive off-diagonal entries as all weights and interpolation coefficients are positive.Yet it is not symmetric: as the vertices used for the interpolation are not incident on the intrinsic edge (i, j) the influence is one-directional.If we want to overcome this asymmetry we need to consider those segments, as in the iARAP discretization.
Note that for symmetry it would be enough to only consider the first segment of an intrinsic edge, not all.That would lead to yet another operator defined by changing the iARAP energy to only depend on the first segment of the intrinsic edge.

Discussion
In this work, we show that the "missing" continuous picture of the original discrete ARAP energy can be found by discretizing the energy of Chao et.al. [CPSS10] over non-conforming elements of dual cells.This interpretation makes clear that the orthogonal dual cell associated with a vertex needs to be immersed.We ensure this by using intrinsic Voronoi cells.The corresponding Laplace-Beltrami operator approximates the intrinsic connectivity using the extrinsic connectivity.It has the same sparsity pattern (stencil) as the standard cotan Laplace-Beltrami operator (i.e. the adjacency matrix of the mesh) but approximates the intrinsic Laplacian.As the connectivity of the intrinsic Laplacian depends on the geometry and may change as the geometry changes, factorizations would have to be recalculated.THe iARAP Laplacian allows for a fixed symbolic factorization or sparse optimizations techniques that assume a fixed sparsity pattern [HTS*22] even after geometry changes.
Our general perspective on the different ARAP discretizations directly leads to further alternatives, which would be worthwhile to investigate in future work.Instead of making the map of the circumcentres part of the optimization, one may assume that they are mapped relative to the vertex positions (i.e.retain their barycentric coordinates).This would lead to another Laplace-Beltrami operator in case of using the extrinsic as well as intrinsic triangulation.
There are various other versions of ARAP modelling, defining the energy in slightly different ways to improve upon the original method.Problems of ARAP also arise from the rigid cells being defined per vertex.The energy penalizes stretching and bending differently because the rotations are defined over vertex stars.Chao et al. [CPSS10] found that the bending penalty decreases as the resolution grows and pointed out that the penalty arises from the vertex wise discretization.This problem can be solved by defining the rotations per triangle and then adding a second term that penalizes bending explicitly [ZG18,LG15].Most of the methods that are based on ARAP or improving on it define energies that sum over edge sets.Those sets contain the edges of a vertex star or a single triangle, whether they define rotations vertex wise or triangle wise.We note that our approach based on the intrinsic Delaunay triangulation is orthogonal to these variations and could be incorporated by using the representation of an intrinsic edge as a sequence of extrinsic segments.It would be interesting to test whether such intrinsic versions of these methods further improve their results on poorly triangulated surfaces.
In a different direction, our derivation also motivates the investigation of other ways to restrict the rotations.Another way to address the above mentioned problems might be a richer function space for the rotations, lifting the restriction from piecewise constant functions.This would also require a proper discretizations of divergence.

Figure 1 :
Figure 1: Deformation of a mesh [SA07, LSHXY22].The constrained vertices (coloured in red) are displaced orthogonal to the plane.The results exhibit the asymmetry in the energy of the spokesand-rims approach.As the initial mesh is Delaunay the iARAP approach coincides with the original ARAP approach.

Figure 2 :
Figure 2: Edges contributing to the ARAP energy associated with vertex i coloured in orange.Line thickness indicates magnitude of the cotan weights -thin lines have weights close to 0.

Figure 3 :
Figure 3: Notation used for primal and dual elements of the mesh.

Figure 4 :
Figure 4: In spokes-and-rims ARAP rotations are fitted per vertex star.The weighting can be interpreted as using only the dark-blue triangles for defining the rotation R i .This means, spokes-and-rims ARAP effectively integrates over 3/4 of the surface.

Figure 5 :
Figure 5: Both the original and spokes-and-rims ARAP assume constant rotations per vertex.The orange indicates the different regions associated to the vertices.Note that the area in the original ARAP is independent of the chosen diagonal, whereas in spokesand-rims ARAP the triangulation is reflected in the cells.

Figure 6 :
Figure6: ARAP[SA07] may be interpreted as the elastic deformation energy of a discontinuous deformation mapping f .The illustration shows how f would affect the input triangulation if it were applied the surface and not only the vertices of the mesh.

Figure 7 :
Figure 7: Notation used for points along the intrinsic polyline.Two adjacent extrinsic triangles where the diagonal edge (u, v) is flipped in the intrinsic triangulation coloured in red.
[JP*18].For the intrinsic Delaunay triangulation we use the signpost data structure [SSC19] implementation of geometry central[SC*19].All results were computed in a single thread on an AMD Ryzen 7 2700X.

Figure 8 :
Figure 8: Visualization of the energy decay per iteration on a mesh with ≈ 20K vertices.The three methods perform roughly similar, with the original ARAP method converging slightly slower.

Figure 9 :
Figure 9: Result for a non-Delaunay mesh.Constraints are coloured in red.The original BUNNY mesh has many ill-shaped triangles in the ears.This results in unwanted deformations for the original ARAP methods.Spokes-and-rims ARAP and iARAP generate similar and better results.

Figure 10 :
Figure 10: The results of spokes-and-rims ARAP (mid) and iARAP (right) are often very similar for reasonably well triangulated surfaces.

Figure 11 :
Figure 11: Deformation of a subdivided cube mesh.The constrained vertices are shifted to the back.The results exhibit the asymmetry in the energy of the spokes-and-rims approach.As the input mesh is Delaunay the result of iARAP and original ARAP coincide.

Figure 12 :
Figure12: A balloon mesh from Thingi10K[ZJ16] is squeezed.Spokes-and-rims ARAP rotates the entire mesh asymmetrically.Original ARAP is omitted as it behaves identical to iARAP on this mesh.

Figure 13 :
Figure13: Deformation of a stool from Thingi10K[ZJ16].Notice the symmetric gap between the legs in the input.As the original ARAP energy diverges we stopped the optimization after a few steps (and the image shows only a subset of the resulting mesh).Spokes-and-rims ARAP and iARAP both converge, but only iARAP preserves the symmetry.

Figure 15 :
Figure 15: Two cases of intrinsic edges intersecting a vertex star coloured in red and orange.Bold arrows follow intrinsic edges.Blue lines are extrinsic edges.Dashed lines are outside of the vertex star of v x .The solid red polyline starts at v x and ends at the boundary of the vertex star of v x .The solid orange polyline starts and ends at the boundary of the vertex star.

©
2023 The Authors.Computer Graphics Forum published by Eurographics -The European Association for Computer Graphics and John Wiley & Sons Ltd.

Table 1 :
Execution times for the precomputation (Cholesky factorization, plus intrinsic Delaunay triangulation for iARAP) and a single iteration (local rotations from extrinsic or intrinsic edges, plus solving) averaged over 100 runs for the different ARAP variants.All times are in milliseconds.

Table 2 :
Execution times for the local step in iARAP depending on the number of segments per intrinsic edge.All instances have the same number of vertices (35K) and edges (105K) but differ in the ratio of segments to edges.Times are in milliseconds.

Table 3 :
Overview of the properties of the two Laplace-Beltrami operators.The operator L is the one used in iARAP, L + is the asymmetric modification.