Controlling Meshes via Curvature: Spin Transformations for Pose-Invariant Shape Processing

TL;DR

This paper introduces discrete spin transformations as a stable, curvature-based method for pose-invariant shape processing of surface meshes, enabling applications like surface fairing and mesh extrusion with minimal distortion.

Contribution

It expands the discrete formulation of spin transformations, allowing controlled curvature manipulation for shape analysis and processing without remeshing, and demonstrates their effectiveness on complex biological structures.

Findings

Stable shape manipulation across large deformations

Effective surface fairing and mesh extrusion

Controlled trade-off between area and angle distortions

Abstract

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.

Equations

$$D_e\varphi =\mu \varphi .$$

$$A_i^{(t)}=A_i(1+h_{i}t+\mathcal{O}(t^2)).$$

$$E_{ij}\triangleq H_{ij}+e_{ij}\in \mathbb{H}.$$

$$E_{ij}$$

$$n_i$$

$$D_{\mathcal{X}}\varphi =\rho A\varphi .$$

$${\tilde{h}}_i{\tilde{A}}_i={\rho }_i{A}_i|{\varphi }_i{|}^2.$$

$$M[q]\triangleq \left[\begin{array}{cccc}[r]a& -b& -c& -d\\ b& a& -d& c\\ c& d& a& -b\\ d& -c& b& a\end{array}\right].$$

$$\underset{\text{upphi}}{\mathrm{argmin}}\underset{D_{\rho }\varphi \simeq 0}{\underbrace{{\text{upphi}}^{\mathsf{T}}({\mathrm{D}}_{\rho }{\mathrm{A}}^{-1}{\mathrm{D}}_{\rho })\text{upphi}}}+\underset{\text{regularization}}{\underbrace{(\text{upphi}-1{)}^{\mathsf{T}}\alpha \mathrm{R}(\text{upphi}-1)}},$$

$$\text{Updelta}\tilde{f}=\nabla \cdot \tilde{e},$$

$$E_{ij}=\frac{|e_{ij}|}{\cos ({\theta }_{ij}/2)}\exp \left(\frac{\pi -{\theta }_{ij}}{2}{u}_{ij}\right).$$

$$L_{\wedge }^2(E)=\mathrm{im}([\nabla \times {]}^{\mathsf{T}})\oplus \underset{\ker (\nabla \times )}{\underbrace{\ker ({\text{Updelta}}_1)\oplus \mathrm{im}(\nabla )}},$$

$$\mathrm{Im}~{}\sum_{i}\bar{\phi}_{i}\Big{(}\sum_{j\sim i}w_{ij}E_{ij}\omega_{ij}^{k}\phi_{j}\Big{)}=0\,,$$

$$2\sum_{i}\dot{\bar{\phi}}_{i}\Big{(}\sum_{j\sim i}w_{ij}E_{ij}\omega_{ij}^{k}\Big{)}=2\sum_{i}A_{i}\dot{\bar{\phi}}_{i}v_{i}^{k}\in\mathbb{R}\,,$$

$$\big{\langle}\delta\!\rho\big{|}z_{\nu}^{k}\big{\rangle}_{0,\mathbb{R}}=0,\quad k=1\cdots b_{1},\quad\nu=\bm{\mathrm{i}},\bm{\mathrm{j}},\bm{\mathrm{k}}\,.$$

$$s_i\triangleq \log A_i/A_i^0-\log ⟨A_i⟩/⟨A_i^0⟩,$$

$${\tilde{E}}_{ij}$$

$$=\overline{(1+{\varphi }_i^{-1}\delta {\varphi }_i)}{E}_{ij}(1+{\varphi }_j^{-1}\delta {\varphi }_j)$$

$$|{\tilde{e}}_{ij}|\cos ({\theta }_{ij}/2)=|1+{\varphi }_i^{-1}\delta {\varphi }_i|\cdot |e_{ij}|\cos ({\tilde{\theta }}_{ij}/2)\cdot |1+{\varphi }_j^{-1}\delta {\varphi }_j|.$$

$${\tilde{A}}_i\approx A_i|1+{\varphi }_i^{-1}\delta {\varphi }_i{|}^4,$$

$${\tilde{s}}_i$$

$$=s_i+\log \frac{|{\varphi }_i+\delta {\varphi }_i{|}^4}{|{\varphi }_i{|}^4}-\log {⟨\frac{|\varphi +\delta \varphi {|}^4}{|\varphi {|}^4}⟩}_A$$

$$\approx s_i+4\left(\frac{⟨{\varphi }_i|\delta {\varphi }_i{⟩}_{\mathbb{H}}}{|{\varphi }_i{|}^2}-{⟨\frac{⟨\varphi |\delta \varphi {⟩}_{\mathbb{H}}}{|\varphi {|}^2}⟩}_A\right)$$

$$\mathrm{F}_{i}=4\Big{(}\langle s\lambda\rangle_{A}\!-\!s_{i}\lambda_{i}\Big{)}\left|\frac{a_{i}}{|\phi_{i}|^{2}}\upphi_{i}\right\rangle_{\mathbb{H}}\,,\quad Q_{i}=\left|\frac{\sqrt{a_{i}\lambda_{i}}}{|\phi_{i}|^{2}}\phi_{i}\right\rangle_{\mathbb{H}}\,,$$

$${\mathrm{L}}_{j1}={|\frac{a_j\sqrt{⟨\lambda {⟩}_A}}{|{\varphi }_j{|}^2}{\varphi }_j⟩}_{\mathbb{H}},{\mathrm{L}}_{j2}={|\frac{a_j{\lambda }_j}{|{\varphi }_j{|}^2}{\varphi }_j⟩}_{\mathbb{H}},{\mathrm{L}}_{j3}={|\frac{a_j}{|{\varphi }_j{|}^2}{\varphi }_j⟩}_{\mathbb{H}}.$$

Full text

1 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 $3$D 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 $\phi\colon\mathcal{F}\rightarrow\mathbb{H}$ (one quaternion per face in the mesh) to a change $\mu\colon\mathcal{F}\rightarrow\mathbb{R}$ in the mean curvature via a first-order differential operator $D_{e}$:

[TABLE]

Typically, a desired change of curvature is specified via $\mu$, yielding a transformation $\phi$ 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 $\mathbb{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.

2 Related work

2.0.1 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).

2.0.2 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.

3 Discrete Geometric Setting

Face edge-constraint nets.

Our work focuses on the case of closed compact orientable surfaces in $\mathbb{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 $\mathcal{G}=(\mathcal{V},\mathcal{F},\mathcal{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)\in\mathcal{E}}$), the normals satisfy the looser condition $n_{i}+n_{j}\perp e_{ij}$. $\mathcal{X}=(\mathcal{G},n)$ is called a face edge-constraint net. For instance standard triangulations with normals orthogonal to faces are face edge-constraint nets.

Discrete mean curvature.

Let $\mathcal{X}$ be a net. $\mathcal{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 $\theta_{ij}$ between planes $P_{i}\!\triangleq\!\mathrm{span}\{n_{i},e_{ij}\}$ and $P_{j}\!\triangleq\!\mathrm{span}\{n_{j},e_{ij}\}$. For standard triangulations, $\theta_{ij}$ is just the bending angle between faces. The integrated mean curvature on edge $e_{ij}$ (Fig. 2(c)) is defined as $H_{ij}\!\triangleq\!|e_{ij}|\tan(\theta_{ij}/2)$. The integrated mean curvature on face $i$ is the sum of its integrated edge curvatures: $H_{i}\!\triangleq\!\sum_{j\in\mathcal{N}(i)}H_{ij}$. The discrete mean curvature $h_{i}\!\triangleq\!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):

[TABLE]

Geometry in the quaternions.

Quaternions $\mathbb{H}$ provide a natural algebraic language for geometry in $\mathbb{R}^{3}$, much like complex numbers for planar geometry. Let $\{1,\bm{\mathrm{i}},\bm{\mathrm{j}},\bm{\mathrm{k}}\}$ denote a basis for $\mathbb{H}$. Elements $v\!\triangleq\!(v_{x},v_{y},v_{z})\!\in\!\mathbb{R}^{3}$ are identified with pure imaginary quaternions $v_{x}\bm{\mathrm{i}}+v_{y}\bm{\mathrm{j}}+v_{z}\bm{\mathrm{k}}\in\mathrm{Im}~{}\mathbb{H}\triangleq\mathrm{span}\{\bm{\mathrm{i}},\bm{\mathrm{j}},\bm{\mathrm{k}}\}$, so that surfaces are naturally immersed in $\mathrm{Im}~{}\mathbb{H}$. Denote by $\bar{q}\triangleq a-(b\bm{\mathrm{i}}\!+\!c\bm{\mathrm{j}}\!+\!d\bm{\mathrm{k}})$ the quaternionic conjugate of $q\!\triangleq\!a\!+\!b\bm{\mathrm{i}}\!+\!c\bm{\mathrm{j}}\!+\!d\bm{\mathrm{k}}\in\mathbb{H}$. The norm $|q|$ of $q$ is defined as the square root of $\bar{q}q=a^{2}\!+\!b^{2}\!+\!c^{2}\!+\!d^{2}$. All $q\neq 0$ admit an inverse $q^{-1}=\bar{q}/|q|^{2}$. Again like complex numbers, quaternions admit a polar decomposition $q\!=\!se^{\theta u}\!=\!s(\cos(\theta)+\sin(\theta)u)$ with $u\!\in\!\mathrm{Im}~{}\mathbb{H}$ a unit vector, which makes their geometric meaning more explicit. Indeed, $v\mapsto qvq^{-1}$, also known as conjugation by $q$, expresses rotation around $u$ by an angle $2\theta$. In the same vein, the expression $\tilde{v}=qv\bar{q}$ conveniently expresses a similarity transformation: $\tilde{v}$ corresponds to $v$ rotated around $u$ by $2\theta$ and rescaled by $s^{2}$.

Hyperedges.

Every edge in the net $\mathcal{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}\in\mathrm{Im}~{}\mathbb{H}$ of the edge:

[TABLE]

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 0.A). With that, spin transformations are introduced in a straightforward manner.

4 Discrete Spin Transformations

Discrete spin transformations.

A discrete spin transformation $\phi$ associates a single quaternion $\phi_{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:

[TABLE]

The elegance of the construct lies in the fact that Eq. (4) does transform a face edge-constraint 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 0.A) and from the constraint on face normals. Furthermore, discrete spin transformations $E\rightarrow_{\phi}\tilde{E}$ are trivially invertible: $\tilde{E}\rightarrow_{\phi^{-1}}\!E$. The integrability condition that each face in the new net closes ($\sum_{j}\!\tilde{E}_{ij}\in\!\mathbb{R}$) is equivalent to the existence of a real valued function $\rho:i\mapsto\rho_{i}\in\mathbb{R}$ over faces such that:

[TABLE]

Equation (5) is the cornerstone of the framework. $D_{\mathcal{X}}$ is henceforth referred to as the intrinsic Dirac operator. $D_{\mathcal{X}}$ sends a quaternionic function over faces to another one such that $(D_{\mathcal{X}}\phi)_{i}\triangleq\sum_{j}E_{ij}\phi_{j}$. Left multiplying both sides by $\bar{\phi}_{i}$, the closedness constraint on faces is immediately apparent: $\bar{\phi}_{i}(D_{\mathcal{X}}\phi)_{i}=\sum_{j}\tilde{E}_{ij}$ must be real-valued. For ease of exposition, the expression in the introduction is formulated using a slightly different yet immediately related operator, the extrinsic Dirac operator $(D_{e}\phi)_{i}\!\triangleq\!\sum_{j}E_{ij}(\phi_{j}-\phi_{i})=(D_{\mathcal{X}}\phi)_{i}-H_{i}\phi_{i}$. 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_{\mathcal{X}}$ creates an explicit relationship between a spin transformation $\phi$ and the discrete (resp. integrated) mean curvature $\tilde{h}$ (resp. $\tilde{H}$) of the new net, namely $\bar{\phi}_{i}(D_{\mathcal{X}}\phi)_{i}=\tilde{H}_{i}\triangleq\tilde{h}_{i}\tilde{A}_{i}$ as long as the new net closes. Coupling with Eq. (5),

[TABLE]

When $\phi\!\coloneqq\!1$ is the identity transform, $\rho_{i}\!=\!h_{i}\!=\!\tilde{h}_{i}$. For smooth $|\phi_{i}|$ and from Eq. (4), $\rho_{i}\sqrt{A_{i}}\!\approx\!\tilde{h}_{i}\sqrt{\tilde{A}_{i}}$. In other words, $\rho_{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 $\mu$ describes a change in half-density instead: $\tilde{h}_{i}\tilde{A}_{i}=(h_{i}+\mu_{i})A_{i}|\phi_{i}|^{2}$.

Dirac operators.

Dirac operators $D_{\mathcal{X}}$ and $D_{e}$ have a number of properties that make them appealing for various tasks in shape analysis. $D_{\mathcal{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 $\mathbb{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 $\mathbb{R}^{3}$ (new transformed $\tilde{\mathcal{X}}$). The first (null) eigenvector of $D_{e}$ is trivial. $D_{\mathcal{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_{\mathcal{X}}$ provides a generally non-trivial immersion with higher smoothness than the original shape (lower Willmore energy $\int|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_{\mathcal{X}}$ (Fig. 3). In this work however, we investigate a strategy closely related to [6] with a fine-grained control over the surface deformations.

5 Algorithms

Remark 5.1

Quaternions $q$ admit representations $M[q]$ as $4\!\times\!4$ real matrices (Eq. (7)), so that standard linear algebra libraries can be used to solve quaternionic linear systems. In particular, $M[\bar{q}]=M[q]^{\mkern-1.5mu\mathsf{T}}$, thus Hermitian (quaternionic) forms are represented by real symmetric matrices. We denote real vectors and matrix representations below with upright bold symbols.

[TABLE]

5.0.1 Overview.

The scalar function $\rho$ 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 $\rho$ can be associated with some $\phi$ such that the integrability condition Eq. (5) is satisfied. Namely, $D_{\rho}\!\triangleq\!D\!-\!\rho$ should have a null eigenvalue. This leads Crane et al. [5] to solve for the smallest eigenvalue $\gamma$ and eigenvector $\phi$, yielding a solution of Eq. (5) up to a small constant shift: $D\phi\!=\!(\rho\!+\!\gamma)\phi$. We propose instead to formulate the objective $D_{\!\rho}\phi\!\simeq\!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 $\phi$ amounts to solving a quadratic problem:

[TABLE]

under a set of linear constraints on $\upphi$. $\bm{\mathrm{A}}$ is a diagonal matrix of face areas. In practice we set $\bm{\mathrm{R}}$ to $\bm{\mathrm{A}}\!+\!\beta\bm{\mathrm{L}}_{f}$, where $\bm{\mathrm{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 $\beta\!\coloneqq\!0$.

Overall, the procedure is as follows: prescribe a scalar function $\rho$ for a target shape or curvature change (sec. 6); then solve for the spin transformation $\phi$ (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.

5.0.2 Computing the new immersion.

Let transformed edges $\tilde{e}_{ij}\!=\!\mathrm{Im}\,\tilde{E}_{ij}$ be indexed by their start and end vertices $v\rightarrow v^{\prime}$. Vertex coordinates $\tilde{f}\colon v\!\in\!V\mapsto\tilde{f}_{v}$ satisfy $\tilde{f}_{v^{\prime}}-\tilde{f}_{v}=\tilde{e}_{v\rightarrow v^{\prime}}$. In practice, we solve the mathematically equivalent (Appendix 0.B) linear system

[TABLE]

where $\Updelta$ and $\nabla\cdot$ 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 ($\tilde{e}$) mesh metric (with empirically identical results). A benefit of working from a discrete setting is that no discretization error is introduced from $\phi$ to the corresponding $\tilde{f}$.

5.0.3 Geometrically constrained flows.

The proposed formulation (Eq. (8)) enables fine-grained control over the flow by prescribing additional constraints. For instance, the method extends to topologies beyond spherical by adding an exactness constraint (Appendix 0.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 $|\phi_{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 $|\phi_{i}|^{4}$ of the spin transformation relates to the local change of area. Thus scale changes $\log{\tilde{A}_{i}/A_{i}}$ (up to global rescaling) can be penalized via a linearized soft constraint over $\phi$ (Appendix 0.D).

5.0.4 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 $\rho$. 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 $\rho$ onto a linear subspace. In particular, Appendix 0.C derives alternative integrability conditions in the form of simple linear constraints on $\rho$, for the proposed discrete geometric framework.

6 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 $20$k). 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).

6.0.1 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)\!\triangleq\!\int_{S}|\nabla f|^{2}dA$ or the Willmore functional $E_{W}(f)\!\triangleq\!\int_{S}h^{2}dA$. Recalling that $\Delta f=h\vec{n}$ and ignoring the dependence of $\Delta$ on $f$, gradient descent on $E_{M}$ (resp. $E_{W}$) yields $\dot{f}\!\propto\!\Delta f$ ($\dot{f}\!\propto\!\Delta^{2}\!f$). The former yields the widespread mean curvature flow $\dot{f}\!\propto\!-h\vec{n}$ that iteratively evolves points along the surface normal $\vec{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 $\dot{h}\!\coloneqq\!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 $\delta h\!\coloneqq\!-\tau h$ is then (optionally filtered and) integrated into a new surface immersion $\tilde{f}$, by computing the corresponding spin transformation as per section 5. Specifically, for a given target curvature $\bar{h}_{i}$ (say $h_{i}\!+\!\delta h_{i}$) and area $\bar{A}_{i}$, we let $\rho_{i}\!\coloneqq\!\bar{h}_{i}\sqrt{\bar{A}_{i}/A_{i}}$ (section 4). The standard unconstrained optimization (Eq. (8)) regularized with the face Laplacian $\bm{\mathrm{L}}_{f}$ (or one of its powers) yields quasi-conformal transformations (Fig. 1). Large steps $\tau\!=\!0.5$–$1$ typically remain stable. Whether $\phi$ is numerically integrable can be checked by monitoring the discrepancy between edges $\tilde{E}$ integrated as per Eq. (4), and edges recomputed from $\tilde{f}$ (after getting $\tilde{f}$ 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 0.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).

6.0.2 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 $\langle h\rangle\vec{n}$, where $\langle h\rangle\!\!=\!\!\int_{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 $\epsilon_{s}=\log{\tilde{A}/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$ to $\tilde{f}$. For a conformal deformation, $Q$ is identically $1$ throughout the mesh. However the area distortion $\epsilon_{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\pm 165$, $\epsilon_{s}\!=\!4.1\pm 1.5$; unconstrained spin transformation: $Q\!=\!1.42\pm 0.08$, $\epsilon_{s}\!=\!2.9\pm 0.3$; area constrained: $Q\!=\!1.7\pm 0.2$, $\epsilon_{s}\!=\!0.85\pm 0.05$. For the area constrained spin transform, the maximum area discrepancy simply reflects a user-specified soft target.

6.0.3 Mesh Extrusion.

The task is now to reconstruct (“extrude”) a shape of interest back from a reference sphere, given its mean curvature $h^{\star}$ and area $A^{\star}$ 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 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 get $\bar{h}_{i}$ from $\delta h_{i}\!\coloneqq\!h_{i}^{\star}\!-\!h_{i}$, and set $\rho_{i}\!\coloneqq\!\bar{h}_{i}\sqrt{A_{i}^{\star}/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.4$mm. 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 ($\sim 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\pm 0.3$mm ($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.

7 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.

Acknowledgments

This work is supported by the EPSRC (grant ref no. EP/P023509/1) and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 757173, project MIRA, ERC-2017-STG). DC is also supported by CAPES, Ministry of Education, Brazil (BEX 1500/15-05). KK is supported by the President’s PhD Scholarship of Imperial College London. IW is supported by the Natural Environment Research Council (NERC).

Appendix 0.A Geometric interpretation of hyperedges

Remark 1

Letting $u_{ij}\!\triangleq\!e_{ij}/|e_{ij}|$ and after straightforward manipulations, we get:

[TABLE]

Thus conjugation by $E_{ij}^{-1}$ sends any vector lying on face $i$ to face $j$, and $-n_{i}$ to $n_{j}$. Intuitively speaking, $E_{ij}$ carries information about a connection structure between the affine spaces of faces $j$ and $i$ [12].

Remark 2

$\bar{E}_{ij}=H_{ij}-e_{ij}=H_{ji}+e_{ji}=E_{ji}$. Moreover, assuming the mesh to be closed, edges sum to [math] over any given face, so that $\sum_{j}E_{ij}=H_{i}\in\mathbb{R}$.

Remark 3

Definition of hyperedges as per Eq. (3) may seem somewhat arbitrary. In fact, $E_{ij}$ is necessarily of the form $|e_{ij}|\tan{\alpha_{ij}}+e_{ij}$, up to a multiplicative constant, under the following mild conditions: (a) the imaginary part of $E_{ij}$ is along $e_{ij}$; (b) $\bar{E}_{ij}=E_{ji}$; and (c) $\sum_{j}\!E_{ij}\!\in\!\mathbb{R}$ iff face $i$ closes. Relating $\alpha_{ij}$ to the bending angle $\theta_{ij}$ is sufficient to guarantee that spin transformations transform a net into another valid net (i.e. the real part $\tilde{H}_{ij}$ of the transformed edge $\tilde{E}_{ij}$ is provably consistent with the constructive definition above).

Appendix 0.B Edge Integration

Let $\nabla\tilde{f}\colon\epsilon\!\triangleq\!(v\!\rightarrow\!v^{\prime})\mapsto(\tilde{f}_{v^{\prime}}\!-\!\tilde{f}_{v})\!\triangleq\!\nabla\tilde{f}(\epsilon)$ the discrete gradient. We are looking for $\tilde{f}$ s.t. $\nabla\tilde{f}(\epsilon)=\tilde{e}_{\epsilon}$. When such an $\tilde{f}$ exists, $\tilde{e}$ is said to be exact (as a discrete $1$-form, defined over edges $\epsilon$). In that case, Eq. (9) follows by taking $\nabla\cdot$ on both sides.

Define the discrete curl operator $[\nabla\!\times\!\tilde{e}](i)\!\triangleq\!\sum_{j\in\mathcal{N}(i)}\tilde{e}_{ij}$, indexing as in section 3. The curl sends a 1-form (over edges) to a 2-form (over faces). If $[\nabla\!\times\!\tilde{e}]$ vanishes everywhere, $\tilde{e}$ is said to be closed. It is easy to verify that $\nabla\!\times\!\nabla\tilde{f}$ is everywhere zero, so that exactness always implies closedness. For a (discrete) spherical topology, the converse holds: closedness implies exactness. Now let $\tilde{e}\!=\!\mathrm{Im}\,\tilde{E}$ be generated by a spin transformation $\phi$ acting on hyperedges $E$ with $e$ closed. Then $\tilde{e}$ is closed iff Eq. (5) is satisfied (immediate from the definition of $D_{\mathcal{X}}$, cf. section 4).

Remark 4

High-level elements of constructive proof are derived from the mesh being simply connected. It is path connected, so we can fix a vertex $v$ and reach any vertex $v^{\prime}$ from $v$ by following a path $\gamma({v\!\rightarrow\!v^{\prime}})$ on edges. Let $\tilde{f}_{v^{\prime}}\!\coloneqq\!\tilde{f}_{v}+\sum_{\epsilon\in\gamma(v\!\rightarrow\!v^{\prime})}\tilde{e}_{\epsilon}$ obtained by summing edges along the path. $\tilde{f}$ is well defined because the value at $v^{\prime}$ is independent of the path. Indeed let $\gamma_{1}$, $\gamma_{2}$ two paths from $v$ to $v^{\prime}$. Following $\gamma_{1}$ then the reverse of $\gamma_{2}$, we run a closed loop. Self-intersections are removed without loss of generality. One can prove by induction on the loop length that edges sum to [math] over the loop, thus the sum over $\gamma_{1}$ and $\gamma_{2}$ are equal. This holds for vertices on a single face by closedness. Closed loops of arbitrary length can always be incrementally shrunk down to this case (by simple connectivity), without changing the sum of edge values (by closedness).

0.B.0.1 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 $\tilde{e}$ can fail to be exact if it has a non-zero harmonic component, i.e. if it can be written as $\tilde{e}\!=\!\omega+\nabla\tilde{f}$ for some discrete 0-form $\tilde{f}\colon V\!\rightarrow\!\mathrm{Im}\,\mathbb{H}$ and harmonic 1-form $\omega\colon E\!\rightarrow\!\mathrm{Im}\,\mathbb{H}$ (s.t. $\omega\!\neq\!0$ is closed and $\Updelta_{1}\omega\!=\!0$). Equivalently $\omega$ is closed with vanishing divergence $\nabla\cdot\omega\!=\!0$. While Eq. (9) still admits a solution $\tilde{f}$, the corresponding edges, $\tilde{e}\!-\!\omega$, are not the prescribed ones. Fortunately, convenient exactness constraints can be derived via the following theorem.

Theorem 0.B.1 (Helmholtz–Hodge decomposition)

The $L^{2}$ space of (alternating) 1-forms (edge flows) $L_{\wedge}^{2}(E)$ on $\mathcal{G}\!\triangleq\!(\mathcal{V},\mathcal{F},\mathcal{E})$ admits an orthogonal decomposition into subspaces of co-exact, harmonic, and exact forms:

[TABLE]

where $\Updelta_{1}\!\triangleq\!-\nabla[\nabla\cdot]+[\nabla\!\times]^{\mkern-1.5mu\mathsf{T}}[\nabla\!\times]$ is the so-called graph Helmholtzian.

Proof

This is Hodge theorem in linear algebra, noting that $\nabla^{\mkern-1.5mu\mathsf{T}}=-[\nabla\cdot]$ and $\nabla\!\times\!\nabla\!=\!0$. Moreover, by Hodge isomorphism theorem, the dimension of $\mathrm{ker}(\Updelta_{1})$ is the $1$st Betty number $b_{1}$, i.e. the number of handles for $\mathcal{G}$ (generally small).

Concretely, let $\omega^{1}\cdots\omega^{b_{1}}$ a set of null eigenvectors for the Helmholtzian $\Updelta_{1}$. Edges $\tilde{e}$ are integrable iff $\tilde{e}$ is closed (Eq. 5) and orthogonal to $\omega^{k}\nu$ ($k\!=\!1\cdots b_{1}$, $\nu\!=\!\bm{\mathrm{i}},\bm{\mathrm{j}},\bm{\mathrm{k}}$) w.r.t. the inner product on $L_{\wedge}^{2}(E)$: $\langle\tilde{e}|\omega^{k}\nu\rangle_{1,\mathbb{H}}\!=\!0$. In Appendix 0.C, these constraints on transformed edges are turned into constraints on the spin transformation $\phi$, or alternatively on the prescribed curvature map $\rho$.

Remark 5 (Cotangent discretization)

Let $\mathcal{X}$ a triangulated net. For $\mathbb{R}$-edge flows $g,\tilde{g}\in L_{\wedge}^{2}(E)$, define the inner product $\langle g|\tilde{g}\rangle_{1,\mathbb{R}}\!\triangleq\!\frac{1}{2}\sum_{\epsilon}w_{\epsilon}g(\epsilon)\tilde{g}(\epsilon)$. Set edge weights $w_{\epsilon}$ to the symmetric expression $w_{v\rightarrow v^{\prime}}\!\triangleq\!\frac{1}{2}(\cot{\angle vv_{1}v^{\prime}}+\cot{\angle v^{\prime}v_{2}v})$ where $v_{1}$ (resp. $v_{2}$) complete the two triangular faces adjacent to the edge $v\!\rightarrow\!v^{\prime}$. On [math]-form, define the standard inner product $\langle f|\tilde{f}\rangle_{0}\!\triangleq\!\sum_{v}f(v)\tilde{f}(v)/A_{v}$ where vertices are weighted by cell areas $A_{v}$. Define the divergence operator $\nabla\cdot\!\triangleq\!-\nabla^{\mkern-1.5mu\mathsf{T}}$ as the negative adjoint of the gradient $\nabla$. Then $\nabla\cdot$ is exactly the cotangent-weighted divergence (the sum of outbound edge flows at $v$ weighted by $w_{v\rightarrow v^{\prime}}$) and $\Updelta=\nabla\cdot\nabla$ the cotangent Laplacian. With this, the relevant inner product $\langle\cdot|\cdot\rangle_{1}$ compatible with the cotangent scheme is now specified.

Appendix 0.C Exactness Constraint

From Appendix 0.B, $E\!\rightarrow_{\phi}\!\tilde{E}$ is integrable if $\tilde{E}$ is closed and orthogonal to $\omega_{k}\nu$, where $\omega^{k}$ spans real-valued harmonic $1$-forms ($k\!=\!1\cdots b_{1}$) and $\nu\!=\!\bm{\mathrm{i}},\bm{\mathrm{j}},\bm{\mathrm{k}}$, w.r.t. the inner product on $L_{\wedge}^{2}(E)$, say $\langle\cdot|\cdot\rangle_{w_{\epsilon}}$ with the notations of Remark 5.

0.C.0.1 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|\mathcal{F}|$ real-valued constraints ($3$ imaginary dimensions, $|\mathcal{F}|$ faces), e.g. by linearizing $\mathrm{Im}(\bar{\phi}_{i}D\phi_{i})\!=\!0$ around the current solution.

0.C.0.2 Exactness.

Exactness is guaranteed if $\mathrm{Re}(\nu\sum_{\epsilon}w_{\epsilon}\tilde{E}_{\epsilon}\omega_{\epsilon}^{k})\!=\!0$. $\omega_{\epsilon}$ and $\mathrm{Im}(\tilde{E})$ 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:

[TABLE]

that can be linearized around the current solution $\phi$. Alternatively, let us derive the corresponding constraint on $\rho$. Consider a time flow $\phi_{t}$, $\rho_{t}$ starting at $\phi_{0}=1$, with time derivative $\dot{\phi}$, $\dot{\rho}$ at $t\!=\!0$. Deriving w.r.t. time, Eq. (5) becomes $D_{e}\dot{\phi}=\dot{\rho}$ and Eq. (12) rewrites as:

[TABLE]

where we use the alternating property to collapse the two sums and define $v^{k}\!\triangleq\!\nabla\!\times\!(wE\omega^{k})$. In other words, for $\nu$ spanning $\mathrm{Im}\,\mathbb{H}$, \big{\langle}\dot{{\phi}}\big{|}v^{k}\nu\big{\rangle}_{0,\mathbb{H}}\!=\!0. Let $z^{k}$ the unique solution to $D_{e}z^{k}\!=\!v^{k}$ and note that $D_{e}(z^{k}\nu)\!=\!(D_{e}z^{k})\nu$. Finally, since $D_{e}$ is self adjoint, \big{\langle}\dot{\phi}\big{|}D_{e}z^{k}\nu\big{\rangle}_{0,\mathbb{H}}=\big{\langle}D_{e}\dot{\phi}\big{|}z^{k}\nu\big{\rangle}_{0,\mathbb{H}} and we obtain:

[TABLE]

where $z_{\nu}^{k}$ are the three imaginary components of $z^{k}$. The constraint can be enforced by projection of the update $\delta\!\rho$ on the orthogonal subspace of the $z_{\nu}^{k}$.

0.C.0.3 To summarize:

1. (i)

Compute the null eigenspace $\omega_{1}\cdots\omega_{b_{1}}$ of the Helmholtzian $\Updelta_{1}\!\triangleq\!-\nabla[\nabla\cdot]+[\nabla\!\times]^{\mkern-1.5mu\mathsf{T}}[\nabla\!\times]$ 2. (ii)

Compute $v^{k}\!\triangleq\!\nabla\!\times\!(wE\omega^{k})$ and $z^{k}$ s.t. $D_{e}z^{k}\!=\!v^{k}$ 3. (iii)

Project $\delta\!\rho$ onto the orthogonal subspace of the imaginary components of $z^{k}$

Appendix 0.D Area distortion

0.D.0.1 Overview.

We wish to penalize local scale changes $\log{A_{i}/A_{i}^{0}}$ (i.e. $A_{i}$ moving away from the initial area distribution $A_{i}^{0}$), relative to the global rescaling $\sum_{i}A_{i}/\sum_{i}A_{i}^{0}$. Therefore the local scale change (with global rescaling factored out) writes as

[TABLE]

with $\langle A_{i}\rangle$ the average face area. We implement a soft constraint of the type $s_{i}^{2}\leq\epsilon^{2}$, where $\epsilon$ defines a tolerance for area distortion. $\epsilon$ can be set by the user or jointly adjusted over the course of the iterations. After introducing Lagrange multipliers $\lambda_{i}$, we are looking at penalties of the form $\sum_{i}\lambda_{i}s_{i}^{2}/2$, which we approximate by linearizing $s_{i}$ w.r.t. a variation $\delta\phi$ of the spin transformation $\phi$. 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.

0.D.0.2 Linearization of $s_{i}$.

We look for a linearized approximation of $\tilde{s}_{i}$ for a change $\delta\phi$ around the spin transformation $\phi$. We start by linearizing Eq. (4). Noting that $\phi_{i}+\delta\phi_{i}=\phi_{i}(1+\phi_{i}^{-1}\delta\phi_{i})$, we get:

[TABLE]

Recalling from Appendix 0.A that $|E_{ij}|\cos{(\theta_{ij}/2)}=|e_{ij}|$ and taking the norm on both sides, we get:

[TABLE]

We can ignore the change in the cosinus of the dihedral angle (to the first order) for simplicity. Secondly we assume $1+\phi^{-1}\delta\phi$ 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 $\delta\phi$ from the change in edge length, and we get:

[TABLE]

where $A_{i}$ is the area when applying $\phi$ to the initial face-edge constraint net, resp. $\tilde{A}_{i}$ when applying $\phi+\delta\phi$. This yields the following expression for $\tilde{s}_{i}$:

[TABLE]

where $\langle\cdot\rangle_{A}$ denotes the spatial average weighted by the face areas $A$, whereas $\langle\cdot|\cdot\rangle_{\mathbb{H}}$ is the inner product on quaternions. In the last expression, we made use of $|\delta\phi_{i}|\ll|\phi_{i}|$, keeping only first order terms.

0.D.0.3 Penalty matrix assembly.

Writing the penalty as $\frac{1}{2}\upphi^{\mkern-1.5mu\mathsf{T}}\bm{\mathrm{Q}}\upphi-\mathrm{F}^{\mkern-1.5mu\mathsf{T}}\upphi$, the penalty matrix $\bm{\mathrm{Q}}$ is the sum of a sparse block diagonal term and $3$ rank-1 terms, $4^{2}(\text{diag}(Q_{i}Q_{i}^{\mkern-1.5mu\mathsf{T}})+\mathrm{L}_{1}\mathrm{L}_{1}^{\mkern-1.5mu\mathsf{T}}-\mathrm{L}_{2}\mathrm{L}_{3}^{\mkern-1.5mu\mathsf{T}}-\mathrm{L}_{3}\mathrm{L}_{2}^{\mkern-1.5mu\mathsf{T}})$. The derivations are tedious but straightforward, yielding:

[TABLE]

[TABLE]

with the use of bra-ket notation, and where $a$ stands for a normalised area $a\!\triangleq\!A/A_{\text{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.