{'fraction_non_alphanumeric': 0.042916615648340886, 'fraction_numerical': 0.008693522713358638, 'mean_word_length': 4.486899563318778, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Controlling Meshes via Curvature: Spin Transformations for Pose-Invariant Shape Processing

We investigate discrete spin transformations, a geometric framework to manipulate surface meshes by controlling mean curvature. Applications include surface fairing -- flowing a mesh onto say, a reference sphere -- and mesh extrusion -- e.g., rebuilding a complex shape from a reference sphere and curvature specification. Because they operate in curvature space, these operations can be conducted very stably across large deformations with no need for remeshing. Spin transformations add to the algorithmic toolbox for pose-invariant shape analysis. Mathematically speaking, mean curvature is a shape invariant and in general fully characterizes closed shapes (together with the metric). Computationally speaking, spin transformations make that relationship explicit. Our work expands on a discrete formulation of spin transformations. Like their smooth counterpart, discrete spin transformations are naturally close to conformal (angle-preserving). This quasi-conformality can nevertheless be relaxed to satisfy the desired trade-off between area distortion and angle preservation. We derive such constraints and propose a formulation in which they can be efficiently incorporated. The approach is showcased on subcortical structures.


Introduction
Generative shape models are tremendously useful in computational anatomy (shape representation, population analysis), medical imaging and computer vision (segmentation, tracking), computer graphics and beyond. Most approaches to statistical shape analysis fundamentally rely on registration, from landmark based representations and active shape models [4,1,25], to medial representations [14] and Principal Geodesic Analysis [9], to deformable registration and diffeomorphometry [8,33]. Registration is known to be a source of bias in shape analysis, but is often a necessary 'evil' because input data does not come pre-aligned in a common reference frame (or pose). In contrast, the shape information of interest is often invariant to the object pose. Our main motivation is to investigate geometric tools that can open the way to learned, pose-invariant generative shape models (specifically, curves and 3D surfaces). The key insight is that mean curvature is pose-invariant and generally characterizes the shape losslessly. This work investigates spin transformations as the algorithmic tool to computationally implement this insight. The cornerstone of the framework lies in a gracefully simple equation that relates a spin transformation φ : F → H (one quaternion per face in the mesh) to a change µ : F → R in the mean curvature via a first-order differential operator D e : Typically, a desired change of curvature is specified via µ, yielding a transformation φ from which a new shape can be constructed. The present work demonstrates this concept and shows its applicability to manipulate (flow and extrude) closed shapes in a stable manner across large deformations. Section 3 reviews the discrete geometric setting, i.e. (i) the geometric objects to which the framework applies, (ii) discrete mean curvature, (iii) background on quaternions as similarity transformations in R 3 . Section 4 introduces discrete spin transformations. Within the framework of spin transformations, the task of flowing a mesh onto a reference shape and that of extruding a shape back from the reference are highly symmetric: both rely on the ability to compute a transformation based on prescribed curvature and area changes. Section 5 gives an overview of the proposed procedure. Section 6 discusses applications and results.

Related work
Pose-invariant shape analysis. Spectral shape descriptors [29,28], built from the spectrum and eigenfunctions of the Laplace(-Beltrami) operator, have achieved popularity in this context, spanning a variety of applications e.g., object retrieval [2], shape dissimilarity quantification [17], analysis of anatomical structures [26,10,31], transfer of structural and functional data [27,20]. Spectral representations pose two challenges: firstly, going back from the spectral descriptor to the corresponding shape is difficult; and secondly, they tend to discard fine-grained, local information in favor of global shape properties and symmetries. To supplement the intrinsic Laplace-Beltrami operator, the extrinsic Dirac operator [19], which carries more information about the shape immersion, has recently been investigated for shape analysis. Geometric deep learning [3] provides the toolset to analyze functions over a fixed graph. It remains unclear how to analyze graphs themselves. The present work contributes with a lossless and invertible mechanism for turning a mesh into a function (the curvature) over a reference template (say, a sphere).
Shape flows, large deformations, conformal maps. Mean curvature flow is the archetypal algorithm for fairing, in part due to its simplicity and intuitive appeal. However mesh quality tends to degrade quickly throughout the flow, requiring tedious monitoring and remeshing to reduce artefacts and prevent singularities [16]. Furthermore it is not suitable for mesh extrusion. Conformal maps are often perceived as the gold standard in such contexts, and spin transformations originate from this perspective [5,6]. Several discretized and discrete quasi-conformal frameworks have been proposed (e.g., [22,18]) on top of an incredibly rich body of theoretical work. Conformal maps have found a natural application in the context of brain mapping [13,11] by mapping the cortical surface to a reference domain. Rather than strictly on conformality, our focus here is on a parametrization of large deformations that (1) works from the shape invariant mean curvature (2) allows to efficiently flow between any shape and a reference. It is more generally related in spirit to large diffeomorphic frameworks [30,21] that can flow a shape from a template and (pose-equivariant) vector field. Our work expands on the framework of discrete spin transformations as introduced by Ye et al. [32]. One of the appeals of a discrete framework is to bypass discretization errors by design and to offer a consistent definition of discrete geometric concepts such as curvature. We introduce the framework to the community and contribute (i) with an optimization strategy that gives finer-grained control over deformations; (ii) by deriving constraints within this formulation for integrability on general topologies, and area preservation; (iii) by exploring its potency for mesh extrusion.

Discrete Geometric Setting
Face edge-constraint nets. Our work focuses on the case of closed compact orientable surfaces in R 3 and follows the discrete geometric setting introduced by Ye et al. [32]. Surfaces are discretized as face edge-constraint nets, generic constructs that encompass but are not restricted to standard triangulated meshes. Let G = (V, F, E) denote the net combinatorics, resp. its vertices, faces and edges. Adjacent faces meet along a single edge. Edges are shared by exactly two adjacent faces. Faces can be arbitrary polygons (such as with simplex meshes [7]). In addition, let each face be assigned a unit normal n, such that for any two adjacent faces i and j joined along edge e ij ((i, j) ∈ E), the normals satisfy the looser condition n i + n j ⊥ e ij . X = (G, n) is called a face edgeconstraint net. For instance standard triangulations with normals orthogonal to faces are face edge-constraint nets. Discrete mean curvature. Let X be a net. X is orientable and, without loss of generality, directed edges e ij are traversed in the direction towards which they point when cycling over vertices of face i. This lets us orient the dihedral angle θ ij between planes P i span{n i , e ij } and P j span{n j , e ij }. For standard triangulations, θ ij is just the bending angle between faces. The integrated mean curvature on edge e ij (Fig. 2(c)) is defined as H ij |e ij | tan(θ ij /2). The integrated mean curvature on face i is the sum of its integrated edge curvatures: H i j∈N (i) H ij . The discrete mean curvature h i H i /A i follows by turning H i into a density over the face. With this, the discrete mean curvature satisfies a discrete counterpart to Steiner's formula. Steiner's formula is a characterization of mean curvature that relates it to the relative change of area when offsetting the surface in the normal direction n by a distance t (replace A i by an infinitesimal area element dA in Eq. (2) for the original formula): Geometry in the quaternions. Quaternions H provide a natural algebraic language for geometry in R 3 , much like complex numbers for planar geometry. Let {1, i, j, k} denote a basis for H. Elements v (v x , v y , v z ) ∈ R 3 are identified with pure imaginary quaternions v x i+v y j+v z k ∈ Im H span{i, j, k}, so that surfaces are naturally immersed in Im H. Denote byq a−(bi+cj+dk) the quaternionic conjugate of q a+bi+cj+dk ∈ H. The norm |q| of q is defined as the square root ofqq = a 2 +b 2 +c 2 +d 2 . All q = 0 admit an inverse q −1 =q/|q| 2 . Again like complex numbers, quaternions admit a polar decomposition q = se θu = s(cos(θ) + sin(θ)u) with u ∈ Im H a unit vector, which makes their geometric meaning more explicit. Indeed, v → qvq −1 , also known as conjugation by q, expresses rotation around u by an angle 2θ. In the same vein, the expressionṽ = qvq conveniently expresses a similarity transformation:ṽ corresponds to v rotated around u by 2θ and rescaled by s 2 .
Hyperedges. Every edge in the net X is associated with a quaternion E ij dubbed hyperedge, with real part equal to the integrated mean curvature H ij at the edge, and imaginary part equal to the embedding e ij ∈ Im H of the edge: Hyperedges are the fundamental structure on which discrete spin transformations act.
They summarize all the geometric information that, along with the mesh combinatorics, allows to reconstruct the discrete surface immersion (Appendix A). With that, spin transformations are introduced in a straightforward manner.

Discrete Spin Transformations
Discrete spin transformations. A discrete spin transformation φ associates a single quaternion φ i with each face i of a face edge-constraint net. The transformation acts on hyperedges E ij and face normals n i as follows: The elegance of the construct lies in the fact that Eq. (4) does transform a face edgeconstraint net into another edge-constraint net. This is easily checked (cf. [32]), with the main elements of the proof stemming from the geometric interpretation of hyperedges (Appendix A) and from the constraint on face normals. Furthermore, discrete spin transformations E → φẼ are trivially invertible:Ẽ → φ −1 E. The integrability condition that each face in the new net closes ( jẼ ij ∈ R) is equivalent to the existence of a real valued function ρ : i → ρ i ∈ R over faces such that: Equation (5) is the cornerstone of the framework. D X is henceforth referred to as the intrinsic Dirac operator. D X sends a quaternionic function over faces to another one such that (D X φ) i j E ij φ j . Left multiplying both sides byφ i , the closedness constraint on faces is immediately apparent:φ i (D X φ) i = jẼ ij must be realvalued. For ease of exposition, the expression in the introduction is formulated using a slightly different yet immediately related operator, the extrinsic Dirac operator It also discards the normalization by A as in [32]. The proposed normalization however mirrors more faithfully the smooth counterpart of the present setting (see e.g. [15]).
The intrinsic Dirac operator D X creates an explicit relationship between a spin transformation φ and the discrete (resp. integrated) mean curvatureh (resp.H) of the new net, namelyφ i (D X φ) i =H i h iÃi as long as the new net closes. Coupling with Eq.
When φ := 1 is the identity transform, ρ i = h i =h i . For smooth |φ i | and from Eq. (4), In other words, ρ i jointly describes the mean curvature and length element. This quantity is precisely known as the mean curvature half-density h|df | in the smooth setting, and is generally in one-to-one correspondence with a given shape. Finally, with the extrinsic Dirac operator, the corresponding µ describes a change in half-density instead: Dirac operators. Dirac operators D X and D e have a number of properties that make them appealing for various tasks in shape analysis. D X and D e are self-adjoint operators. Dirac operators relate to square roots of the Laplace-Beltrami operator L. Whereas L captures the intrinsic manifold geometry and is invariant by isometry, the Dirac operators can disambiguate much more about the surface immersion into R 3 . We refer the reader to [19,32] for a discussion from this perspective. The eigenvectors of Dirac operators all satisfy Eq. (5) and thus provide new immersions of the abstract manifold into R 3 (new transformedX ). The first (null) eigenvector of D e is trivial. D X cannot have a null eigenvalue for closed surfaces (e.g. spherical topology) of practical interest in the present work, since that would result in a minimal closed surface with everywhere zero mean curvature. The smallest eigenvector of D X provides a generally non-trivial immersion with higher smoothness than the original shape (lower Willmore energy |h| 2 dA). Ye et al. [32] explore this mechanism for the purpose of surface fairing. The next leading eigenvectors give some geometric insight into D X (Fig. 3). In this work however, we investigate a strategy closely related to [6] with a fine-grained control over the surface deformations.

Algorithms
Remark 51. Quaternions q admit representations M [q] as 4×4 real matrices (Eq. (7)), so that standard linear algebra libraries can be used to solve quaternionic linear systems. In particular, M [q] = M [q] T , thus Hermitian (quaternionic) forms are represented by real symmetric matrices. We denote real vectors and matrix representations below with upright bold symbols.
Overview. The scalar function ρ introduced in section 4 provides the primary degrees of freedom for mesh manipulation, and it tightly relates to mean curvature. Of course only a subset of functions ρ can be associated with some φ such that the integrability condition Eq. (5) is satisfied. Namely, D ρ D −ρ should have a null eigenvalue. This leads Crane et al. [5] to solve for the smallest eigenvalue γ and eigenvector φ, yielding a solution of Eq. (5) up to a small constant shift: Dφ = (ρ+γ)φ. We propose instead to formulate the objective D ρ φ 0 as a minimization problem. This gives fine-grained control to add specifications (e.g. smoothness, area distortion), many of which can be efficiently expressed as linear(ized) constraints or quadratic regularizers, within a unified formulation. Thus finding φ amounts to solving a quadratic problem: under a set of linear constraints on φ. A is a diagonal matrix of face areas. In practice we set R to A+βL f , where L f is an integrated Laplacian over faces. The eigensystem actually solved in [5] closely relates to the simplest case where there are no constraints and β := 0.
Overall, the procedure is as follows: prescribe a scalar function ρ for a target shape or curvature change (sec. 6); then solve for the spin transformation φ (Eq. (8)); finally compute new hyperedges (Eq. (4)) and solve a linear system for the new vertex coordinates (Eq. (9)). The steps are typically iterated over, resulting in a flow.
Computing the new immersion. Let transformed edgesẽ ij = ImẼ ij be indexed by their start and end vertices In practice, we solve the mathematically equivalent (Appendix B) linear system where ∆ and ∇· are the standard discrete (cotangent) mesh Laplacian and divergence operators [23]. This method of integration is robust to numerical errors. The Laplacian and divergence are computed w.r.t. either the source (e) or target (ẽ) mesh metric (with empirically identical results). A benefit of working from a discrete setting is that no discretization error is introduced from φ to the correspondingf .
Geometrically constrained flows. The proposed formulation (Eq. (8)) enables finegrained control over the flow by prescribing additional constraints. For instance, the method extends to topologies beyond spherical by adding an exactness constraint (Appendix C). The mapping can also be encouraged to preserve angles (i.e. conformality) and minimize area distortion. Conformality is key in preserving mesh quality across exceptionally large deformations, which prevents considerable loss of numerical stability. It is intuitively described as circles being locally transformed into circles, or indeed texture-preserving (Fig. 1). Quasi-conformality is inherent to the present framework. From Eq. (4), the relative length of edges is preserved as soon as |φ i | varies smoothly across faces. On the other hand large area distortion can be introduced, particularly in regions of high curvature. In some applications, we may prefer to trade off some distortion of angles for a better preservation of areas. We note again from Eq. (4) that the magnitude |φ i | 4 of the spin transformation relates to the local change of area. Thus scale changes logÃ i /A i (up to global rescaling) can be penalized via a linearized soft constraint over φ (Appendix D). Filtering in curvature space. As described in [6] in a related setting, the flow of the spin transformation can also be altered by directly manipulating ρ. The rate of change for geometric features of various scales can be tweaked by manipulating its frequency spectrum. Moreover some constraints can be efficiently enforced by orthogonal projection of ρ onto a linear subspace. In particular, Appendix C derives alternative integrability conditions in the form of simple linear constraints on ρ, for the proposed discrete geometric framework.

Applications
This section showcases the approach on a collection of structured meshes of subcortical structures from the UK Biobank database [24]. The typical mesh size is of a few thousand nodes (up to 20k). The framework was implemented in numpy. The tool mostly relies on efficient (sparse) linear algebra. Experiments were run on a standard laptop (i7-8550U CPU @ 1.80GHz).
Surface Fairing. Surface fairing is the process of producing successively smoother approximations of a mesh geometry f . Most algorithms proceed by minimizing a fairing energy, such as the membrane energy E M (f ) S |∇f | 2 dA or the Willmore functional E W (f ) S h 2 dA. Recalling that ∆f = hn and ignoring the dependence of ∆ on f , gradient descent on E M (resp. E W ) yieldsḟ ∝ ∆f (ḟ ∝ ∆ 2 f ). The former yields the widespread mean curvature flowḟ ∝ −hn that iteratively evolves points along the surface normal n with a magnitude proportional to the mean curvature h. Crane et al. [6] first suggested in the context of spin transformations to optimize E W directly w.r.t. h, yielding the simple flowḣ := h in curvature space. A benefit of the approach is to decouple time and spatial integration, yielding numerically stable solutions across large time steps. We follow the same strategy. The prescribed change of curvature δh := −τ h is then (optionally filtered and) integrated into a new surface immersionf , by computing the corresponding spin transformation as per section 5. Specifically, for a given target curvatureh i (say h i +δh i ) and areaĀ i , we let ρ i :=h i Ā i /A i (section 4). The standard unconstrained optimization (Eq. (8)) regularized with the face Laplacian L f (or one of its powers) yields quasi-conformal transformations (Fig. 1). Large steps τ = 0.5-1 typically remain stable. Whether φ is numerically integrable can be checked by monitoring the discrepancy between edgesẼ integrated as per Eq. (4), and edges recomputed from f (after gettingf from Eq. (9)). The closedness generally holds within a few percent across several large steps without an explicit constraint, and within 10 −6 with an explicit constraint (Appendix B). A trade-off between conformality and area distortion is achieved by weighing in a soft constraint on the square norm of the logarithmic area distortion (Fig. 4).
Comparison to Mean Curvature Flow. The procedure is compared with an incompressible mean curvature flow. Incompressibility is enforced to make the flow less prone to develop singularities, by adding a balloon energy h n, where h = S hdA is the average mean curvature. Two metrics of interest, defined over the mesh surface, are the conformality error Q and the logarithmic area distortion s = logÃ/A (after normalising to the same total area). The quality factor Q measures how close-to-conformal a transformation is, as the ratio of the largest to smallest eigenvalues of the Jacobian of the mapping from f tof . For a conformal deformation, Q is identically 1 throughout the mesh. However the area distortion s may become significant. Fig. 4 exemplifies the general observation that the mean curvature flow realises a suboptimal trade-off between angle and area preservation. As expected, unconstrained discrete spin transformations are quasi-conformal. Unavoidable area distortion is introduced but mesh elements retain their original quality (right column, top and middle). To contrast, the mean curvature flow arbitrarily destroys the mesh quality, angle and area ratios in regions of high curvature. Area constrained discrete spin transformations implement a sensible compromise, whereby (i) area distortion is lessened; (ii) numerical stability is preserved; (iii) the conformal error increases rather uniformly over the entire mesh, leading to a graceful, slower loss of mesh quality. For surface fairing to a sphere, averaged over a random subset of 100 meshes in the dataset and taking the maximum over the mesh surface, we get the following -mean curvature flow: Q = 97 ± 165, s = 4.1 ± 1.5; unconstrained spin transformation: Q = 1.42 ± 0.08, s = 2.9 ± 0.3; area constrained: Q = 1.7 ± 0.2, s = 0.85 ± 0.05. For the area constrained spin transform, the maximum area discrepancy simply reflects a user-specified soft target. Mesh Extrusion. The task is now to reconstruct ("extrude") a shape of interest back from a reference sphere, given its mean curvature h and area A mapped onto the sphere surface. There is to our knowledge very little done in that direction, even in related work [6,32]. To emphasize, we only wish to recover the original mesh up to pose Fig. 5. Example extrusion of a brain stem from the reference sphere. The original shape is overlaid as a wireframe. (1st and 2nd) Close to the initial stage. The target mean curvature map is displayed, rather than the reconstruction error. Note that the shape flow in the next stages intuitively matches the information captured in these maps. (3rd and 4th) Intermediate stages in the flow, with overlayed reconstruction error. (5th to last) Reconstructed mesh from two views. and scale. Encoding scale presents little difficulty, and shape is invariant under changes of pose. To evaluate the reconstruction accuracy, we rigidly align and rescale the extruded shape to the original one, and compute the maximum distance from points on the extruded mesh to the original surface. The strategy for extrusion closely mirrors that of mesh fairing, whereby we geth i from δh i := h i −h i , and set ρ i :=h i A i /A i . As a preliminary comment, note that the degree of challenge regarding mesh extrusion critically depends on the exact experimental setting and goal, as contrasted in the two following settings. The first experiment aims to estimate the accuracy that can be reached in the best scenario (somewhat upper bounded by the registration error). We take a collection of 300 subcortical meshes from the UK Biobank (incl. brain stems, caudate, putamen, accumbens, amygdala, hippocampus, thalamus, palladium) and flow them onto the unit sphere. We do not perform remeshing, only interpolating relevant maps to nodes and back to faces. We then directly reconstruct the mesh as described above. On average over the dataset, the maximum point-to-surface error is of 0.4mm. The distribution of error is widely spread over different structures, the most challenging being caudates (1.4) and hippocampi (1.2); and the least ones being the accumbens, amygdala, palladium and thalamus (∼ 0.01-0.02). This matches our expectations, given that caudates and hippocampi are in fact highly non spherical. Thus very significant area or angle distortion is introduced when mapping onto the sphere. The second experiment investigates a more challenging setup, whereby the flowed surface is remapped onto a reference sphere with uniform meshing. Shape-specific vertex density as well as face aspect ratio, which reflect the area and angle distortion introduced during the fairing, are thus discarded. We experiment with a set of 100 brain stems (Fig. 5), which represent a happy medium between the most challenging and trivial structures, with a maximum reconstruction error of 1.4 ± 0.3mm (2-4%). For the most challenging structures, various strategies to guide the reconstruction using either additional information obtained during the flow, or multiscale approaches with hierarchical encoding could be considered. This is left to explore in future work.

Conclusion
We have presented a method to manipulate surface meshes across very large deformations by prescribing mean curvature (half-density). The framework is well suited for mesh fairing and extrusion, e.g. to map shapes to, or back from a unit sphere. As a perspective, we believe the approach to have potential for pose-invariant shape analysis, specifically for generative modeling. Indeed mean curvature together with the metric generally is in one-to-one correspondence with the (closed) shape; this is in particular true for a spherical topology. We have shown how spin transformations computationally implement this insight. Therefore the shape geometry could be losslessly encoded as a scalar function on a template, making the modeling task more amenable to learning. In the smooth setting, spin transformations are a subgroup of conformal maps. This partly explains their numerical stability across large flow steps, a property inherited in the discrete setting. However, conformal maps can introduce significant area distortion, e.g. when flowing highly curved objects. An advantage of discrete spin transformations is to relax exact conformality, and allow the user to trade off angle for area preservation.
Non-simply connected topologies. Consider a path connected mesh, but possibly with handles (note that closed loops circling a handle cannot be continuously shrunk down to a trivial loop). A closed 1-formẽ can fail to be exact if it has a non-zero harmonic component, i.e. if it can be written asẽ = ω+∇f for some discrete 0-formf : V → Im H and harmonic 1-form ω : E → Im H (s.t. ω = 0 is closed and ∆ 1 ω = 0). Equivalently ω is closed with vanishing divergence ∇ · ω = 0. While Eq. (9) still admits a solutioñ f , the corresponding edges,ẽ−ω, are not the prescribed ones. Fortunately, convenient exactness constraints can be derived via the following theorem.
Theorem 1 (Helmholtz-Hodge decomposition). The L 2 space of (alternating) 1-forms (edge flows) L 2 ∧ (E) on G (V, F, E) admits an orthogonal decomposition into subspaces of co-exact, harmonic, and exact forms: where is the so-called graph Helmholtzian.
Proof. This is Hodge theorem in linear algebra, noting that ∇ T = −[∇·] and ∇×∇ = 0. Moreover, by Hodge isomorphism theorem, the dimension of ker(∆ 1 ) is the 1st Betty number b 1 , i.e. the number of handles for G (generally small).
Remark 5 (Cotangent discretization). Let X a triangulated net. For R-edge flows g,g ∈ L 2 ∧ (E), define the inner product g|g 1,R 1 2 w g( )g( ). Set edge weights w to the symmetric expression w v→v 1 2 (cot ∠vv 1 v + cot ∠v v 2 v) where v 1 (resp. v 2 ) complete the two triangular faces adjacent to the edge v → v . On 0-form, define the standard inner product f |f 0 v f (v)f (v)/A v where vertices are weighted by cell areas A v . Define the divergence operator ∇· −∇ T as the negative adjoint of the gradient ∇. Then ∇· is exactly the cotangent-weighted divergence (the sum of outbound edge flows at v weighted by w v→v ) and ∆ = ∇ · ∇ the cotangent Laplacian. With this, the relevant inner product ·|· 1 compatible with the cotangent scheme is now specified.
Closedness. This is Eq. (5) and already core to the framework. Algorithmically, Eq. (8) only guarantees the closedness to approximately hold, but we have observed very good agreement in practice without further constraint. If necessary, closedness can be strictly enforced as a set of 3|F| real-valued constraints (3 imaginary dimensions, |F| faces), e.g. by linearizing Im(φ i Dφ i ) = 0 around the current solution.
Exactness. Exactness is guaranteed if Re(ν w Ẽ ω k ) = 0. ω and Im(Ẽ) are alternating (e.g. w ij = −w ji with the conventions of section 3) so this rewrites as a set of b 1 (3b 1 real-valued) constraints: that can be linearized around the current solution φ. Alternatively, let us derive the corresponding constraint on ρ.
where we use the alternating property to collapse the two sums and define v k ∇ × (wEω k ). In other words, for ν spanning Im H, φ v k ν 0,H = 0. Let z k the unique solution to D e z k = v k and note that D e (z k ν) = (D e z k )ν. Finally, since D e is self adjoint, φ D e z k ν 0,H = D eφ z k ν 0,H and we obtain: where z k ν are the three imaginary components of z k . The constraint can be enforced by projection of the update δρ on the orthogonal subspace of the z k ν .

−∇[∇·] + [∇×] T [∇×]
(ii) Compute v k ∇×(wEω k ) and z k s.t. D e z k = v k (iii) Project δρ onto the orthogonal subspace of the imaginary components of z k D Area distortion Overview. We wish to penalize local scale changes log A i /A 0 i (i.e. A i moving away from the initial area distribution A 0 i ), relative to the global rescaling i A i / i A 0 i . Therefore the local scale change (with global rescaling factored out) writes as with A i the average face area. We implement a soft constraint of the type s 2 i ≤ 2 , where defines a tolerance for area distortion. can be set by the user or jointly adjusted over the course of the iterations. After introducing Lagrange multipliers λ i , we are looking at penalties of the form i λ i s 2 i /2, which we approximate by linearizing s i w.r.t. a variation δφ of the spin transformation φ. We have found this mechanism to hold over large integration steps in practice. This can be better intuited by looking at the nature of the approximations made during the linearization (see below). The approximate quadratic energy is the sum of a sparse block diagonal matrix and a low-rank (dense) term. Woodbury matrix identities allow to solve quadratic systems involving this energy without directly storing or manipulating the dense matrix.
Linearization of s i . We look for a linearized approximation ofs i for a change δφ around the spin transformation φ. We start by linearizing Eq. (4). Noting that φ i +δφ i = φ i (1 + φ −1 i δφ i ), we get: Recalling from Appendix A that |E ij | cos (θ ij /2) = |e ij | and taking the norm on both sides, we get: We can ignore the change in the cosinus of the dihedral angle (to the first order) for simplicity. Secondly we assume 1 + φ −1 δφ to be close to conformal. This assumption is coherent with the spirit of the framework, and should hold regardless if mesh quality is to be preserved locally in time. The overall transformation is still expected to progressively drift from quasi-conformality to accommodate area preservation. With this we can approximate the change in area for a small variation δφ from the change in edge length, and we get:Ã where A i is the area when applying φ to the initial face-edge constraint net, resp.Ã i when applying φ + δφ. This yields the following expression fors i : = s i + log where · A denotes the spatial average weighted by the face areas A, whereas ·|· H is the inner product on quaternions. In the last expression, we made use of |δφ i | |φ i |, keeping only first order terms.
Penalty matrix assembly. Writing the penalty as 1 2 φ T Qφ − F T φ, the penalty matrix Q is the sum of a sparse block diagonal term and 3 rank-1 terms, . The derivations are tedious but straightforward, yielding: with the use of bra-ket notation, and where a stands for a normalised area a A/A tot . The sum of rank-1 updates might be degenerate (for instance it is rank-1 if all multipliers are equal). SVD decomposition can be used to derive an equivalent, non degenerate low-rank basis of vectors.