Total Variation of the Normal Vector Field as Shape Prior with Applications in Geometric Inverse Problems

An analogue of the total-variation prior for the normal vector field along the boundary of smooth and non-smooth shapes in 3D is introduced. Its analysis in the smooth case is based on a differential geometric setting in which the unit normal vector is viewed as an element of the two-dimensional sphere manifold. This novel functional is subsequently extended to piecewise flat, triangulated surfaces as they occur for instance in finite element computations. The ensuing discrete functional agrees with the discrete total mean curvature known in discrete differential geometry. A split Bregman iteration is proposed for the solution of shape optimization problems in which the total variation of the normal appears as a regularizer. Both the continuous and discrete settings are detailed. Unlike most other priors such as surface area, the new functional allows for piecewise flat shapes in the discrete setting. As an application, a geometric inverse problem of inclusion detection type is considered. Numerical experiments confirm that polyhedral shapes can be identified quite accurately.


Introduction
The total variation (TV) functional is popular as a regularizer in imaging and inverse problems; see for instance Rudin, Osher, Fatemi, 1992;Chan, Golub, Mulet, 1999;Bachmayr, Burger, 2009;Langer, 2017 andVogel, 2002, Chapter 8.For a real-valued function u ∈ W 1,1 (Ω) on a bounded domain Ω, the total variation seminorm is defined as (1.1) Notice that we restrict the discussion to the isotropic case here, i.e., | • | 2 denotes the Euclidean norm.(1.1) extends to less regular functions whose distributional gradient exists only in the sense of measures.We refer the reader to Giusti, 1984;Attouch, Buttazzo, Michaille, 2006 for an extensive discussion of functions of bounded variation.The utility of (1.1) as a regularizer, or prior, lies in the fact that it favors piecewise constant solutions.
In this paper we introduce a novel regularizer based on the total variation and demonstrate its utility in geometric inverse problems.In the latter class, the unknown which one seeks to recover is a shape Ω ⊂ R 3 , which might represent the location of a source or inclusion inside a given, larger domain, or the geometry of an inclusion or a scatterer.The boundary of Ω will be denoted by Γ.
The novel functional, which we term the total variation of the normal field along Γ, is defined by In (1.2), n is the outer unit normal vector field along Γ, i.e., n belongs to the manifold S 2 = {v ∈ R 3 : |v| 2 = 1} pointwise.Moreover, D Γ n denotes the derivative (push-forward) of the normal vector field, and {ξ 1 (s), ξ 2 (s)} denotes an orthonormal basis (w.r.t. the Euclidean inner product in the embedding Γ ⊂ R 3 ) of the tangent spaces T s Γ along Γ.Finally, | • | g denotes the norm induced by a Riemannian metric on S 2 .We will consider the metric induced from embedding S 2 in R 3 , i. e. the distance induced by this metric is the arc length distance and the curvature is 1.We write | • | g for the norm induced by the Riemannian metric g on the tangent spaces S 2 .
A thorough introduction to |n| T V (Γ) and its properties will be given in Section 2. Nevertheless we wish to point out already at this point a number of properties of (1.2) which set it apart from (1.1): (1) The variable on which (1.2) depends is the domain Ω.Since the normal vector field n in turn depends on Ω, both the integration domain Γ and the integrand in (1.2) depend on Ω.By contrast, Ω is fixed in (1.1), where u is the variable.
(2) The normal vector field, whose pointwise variation the total variation functional (1.2) seeks to capture, is manifold-valued with values in S 2 .By contrast, the function u in (1.1) is real-valued.
(3) It is well known that the TV functionals penalize jumps and non-zero gradients.Consequently, the minimization of (1.1) avoids unnecessary variations of u and thus favors piecewise constant minimizers.The situation is slightly different for (1.2) since we are considering closed surfaces Γ, which yields a periodicity constraint for the normal vector field n.In this setting, unnecessary variations of n correspond to non-convex regions of the enclosed body Ω.Consequently, the minimization of (1.2) favors convex shapes and, more precisely, spheres; see Theorem 2.4.
Further properties of (1.2) and its discrete counterpart will be discussed in Sections 2 and 4, respectively.
In this paper we are interested in utilizing the total variation of the normal (1.2) as a prior in geometric inverse problems involving a partial differential equation (PDE).As an example, we discuss an inclusion detection problem motivated by geophysical applications.The aforementioned and many other problems of interest can be cast in the form Minimize (u(Ω), Ω) + β |n| T V (Γ) w.r.t.Ω in a suitable class of domains. (1.3) Here u(Ω) denotes the solution of the problem specific PDE, which depends on the unknown domain Ω.Moreover, represents a loss function, such as a least squares function.
The coupling between Ω and its normal vector field n makes the minimization of (1.3) algorithmically challenging.Moreover, since the integrand in (1.2) is zero on flat regions (with constant normal) of Γ, (1.2) and thus (1.3) cannot be expected to be shape differentiable, although the first (loss function) part pertaining to the PDE often is.We therefore resort to a splitting approach in the spirit of Goldstein, Osher, 2009, where d = ∇u was introduced as an independent variable in the context of the total variation functional (1.1).The variables u and d are coupled through a constraint, which is then handled in an Alternating Direction Method of Multipliers (ADMM) framework.We refer the reader to Glowinski, Marroco, 1975;Goldstein, Bresson, Osher, 2010;Goldstein, O'Donoghue, et al., 2014 for more on ADMM.
In our setting, the role of u in (1.1) is played by the normal vector field n depending on the shape Ω, and thus ∇u is replaced by derivative D Γ n, also known as the pushforward of n.In our proposed splitting, we introduce a new variable d, independent of Ω and its normal vector field n, and require the coupling condition d = D Γ n to hold across Γ.For piecewise flat surfaces, such as those represented by simplicial meshes utilized in computations, the normal vector n is constant on facets.As will become clear further below, in this case the push-forward D Γ n needs to be replaced by the manifold analogue of its jump, i.e., the logarithmic map in the appropriate tangent space of S 2 .An outstanding feature of the proposed splitting is that the two subproblems, the minimization w.r.t.Ω and w.r.t.d, are directly amenable to numerical algorithms.The former is a smooth shape optimization problem, and the latter turns out to be solvable explicitly as a shrinkage problem in the respective tangent spaces.
Regularizing functionals involving the normal vector for shape and geometry optimization problems have been considered elsewhere in the literature, mainly in the context of mesh denoising and surface fairing.Since some background in differential geometry will be required in order to allow a comparison with our novel functional (1.2) we postpone the discussion of related work to Section 2.
Although many optimization algorithms have been recently generalized to Riemannian manifolds, the split Bregman method for manifolds proposed in this paper is new to the best of our knowledge.For a general overview of optimization on manifolds, we refer the reader to Absil, Mahony, Sepulchre, 2008.Examples of non-smooth methods for manifold-valued total variation problems have been introduced for instance in Lellmann et al., 2013.A discrete graph-based setting has been considered in Bergmann, Tenbrinck, 2018, which includes a discrete model similar to (4.3).A splitting scheme, the so-called half-quadratic minimization, was introduced by Bergmann, Chan, et al., 2016.
The structure of the paper is as follows.In the following section we provide an analysis of (1.2) and its properties.We also compare (1.2) to geometric functionals appearing elsewhere in the literature.Section 3 is devoted to the formulation of an ADMM method which generalizes the split Bregman algorithm to the manifoldvalued problem (1.3).In Section 4 we address a discrete counterpart of (1.2) on simplicial meshes, as they are frequently used in finite element discretizations of PDEs.The discussion of the corresponding discrete ADMM method is given in Section 5. Section 6 is devoted to the description of implementation details for an inclusion detection problem of type (1.3), motivated by geophysical application, in the finite element framework FEniCS.Corresponding numerical results are presented in Section 7.

Analysis of the Total Variation of the Normal
In this section we discuss (1.2) in detail and relate it to other geometric functionals used previously in the literature.A minimal background in differential geometry of surfaces is required, which we recall here and refer the reader to do Carmo, 1976;Gray, Abbena, Salamon, 2006;Kühnel, 2013 for a thorough introduction.
2.1.Preliminaries.Throughout this section we assume that the boundary Γ of the unknown bounded domain Ω is a smooth, orientable manifold of dimension 2 without boundary, embedded in R 3 .Therefore we can think of tangent vectors at s ∈ Γ to be elements of the appropriate two-dimensional subspace (the tangent plane) of R 3 .This tangent plane at s is denoted by T s Γ.Each tangent plane is endowed with the Riemannian metric furnished by the embedding via the pull-back of the Euclidean metric in R 3 .In other words, the inner product of two vectors In what follows, {ξ 1 (s), ξ 2 (s)} denotes an orthonormal basis in T s Γ.As the following lemma shows, the choice of this basis and how it varies with s will not matter.
Outward pointing unit normal vectors n along Γ will be considered elements of the two-dimensional smooth manifold S 2 .The derivative or push-forward of the normal map n is denoted by In what follows, we will suppress the dependence on the point s ∈ Γ where possible.
Lemma 2.1.The total variation of the normal (1.2) is independent of the choice of the orthonormal basis in the tangent spaces T s Γ.
Proof.Consider a point s ∈ Γ and suppose that {ξ 1 , ξ 2 } and {η 1 , η 2 } are two orthogonal bases of T s Γ.Then there exists an orthogonal matrix Then the integrand in (1.2) can be written as g .This concludes the proof.
Similarly as we do for Γ, we consider S 2 embedded into R 3 and therefore we can conceive the tangent space T n(s) S 2 as a two-dimensional plane in R 3 tangent to the sphere S 2 .We endow T n(s) S 2 with the Riemannian metric furnished by the pull-back of the Euclidean metric as well, which we denote by g(•, •) to distinguish it from the Riemannian metric on T s Γ.In fact, T n(s) S 2 is clearly parallel to T s Γ, see Figure 2.1.We can therefore identify the two tangent spaces and we write Let us argue that (1.2) generalizes (1.1).Since the normal vector field n replaces the scalar-valued function u in (1.1), assume for the moment that n maps into R instead of S 2 .Then the tangent space T n R is equal to R, endowed with its usual inner product.Finally, the manifold Γ in (1.2) takes the role of Ω ⊂ R 2 in (1.1).By Lemma 2.1, we can choose, without loss of generality, the basis ξ 1 = (1, 0) and ξ 2 = (0, 1) .Consequently, (1.2) becomes 2.2.Relation to Curvature.In order to relate (1.2) with regularizing geometric functionals appearing elsewhere in the literature, we take a second look at the integrand.To this end, we recall that the normal field operator N Γ : Γ → S 2 is also known as the Gauss map; see for instance Kühnel, 2013, Chapter 3. Its derivative at s ∈ Γ maps tangent directions in T s Γ into tangent directions in T n(s) S 2 ∼ = T s Γ.
With the latter identification, the derivative of the Gauss map is known as the shape operator S : T s Γ → T s Γ.
Theorem 2.2.The integrand in (1.2) satisfies Proof.Consider the square of the integrand, ). Due to Lemma 2.1 we can choose ξ 1 , ξ 2 to be normalized eigenvectors in T s Γ corresponding to the eigenvalues k 1 , k 2 , respectively.Therefore we get 2.3.Comparison with Preliminary Work.The representation of the integrand from Theorem 2.2 allows us to rewrite (1.2) as the integral over the root mean square curvature, and compare it with related functionals appearing in the literature.The quantity is known as the integral over the total curvature (although this term is also used for other quantities in the literature).The functional (2.4) has a long tradition in surface fairing applications and can be interpreted as a surface strain energy, see for instance Lott, Pullin, 1988;Hagen, Schulze, 1987;Welch, Witkin, 1992;Halstead, Kass, DeRose, 1993;Welch, Witkin, 1994;Greiner, 1994;Tasdizen et al., 2003. Since (2.4) corresponds to Ω |∇u| 2 2 dx in imaging applications, which leads to a Laplacian in the associated optimality conditions (and thus also in the corresponding L 2 -gradient flow), (2.4) tends to smooth the surface and its features.
By contrast, the functional (2.3) seems to have made very few appearances in the mathematical literature.We are aware of the PhD thesis Maekawa, 1993, Chapter 6 and the subsequent book publication Patrikalakis, Maekawa, 2001 where it was used to guide mesh generation.In Ateshian, Rosenwasser, Mow, 1992;Marzke et al., 2012 the pointwise root mean square curvature is used as a measure of flatness in biomedical classification problems, in Pulla, Razdan, Farin, 2001 for the purpose of surface segmentation and in Wu, Ma, et al., 2010 it is used as an aid to visualize vascular structures.We also mention that the logarithm of the root mean square curvature is known as the curvedness and it plays a role in the classification of intermolecular interactions in crystals; see for instance McKinnon, Spackman, Mitchell, 2004.We are however not aware of any use of (1.2) or its equivalent form (2.3) as a prior in geometric inverse problems.
We regard (1.2) as a natural extension of the total-variation seminorm (1.1) to measure surface flatness, but other extensions are certainly possible.Notably, the authors in Elsey, Esedoglu, 2009 propose the total absolute Gaussian curvature for the same purpose.From the Gauss-Bonnet theorem (see for instance Gray, Abbena, Salamon, 2006, Chapter 27 or Kühnel, 2013, Chapter 4F) it follows that the boundaries Γ of convex domains Ω will be the global minimizers of (2.5), and they yield a value of 4π.Thus (2.5) promotes domains which are "as convex as possible".
It should be noted that the classical total variation seminorm (1.1) is not invariant with respect to scale.In fact, it is easy to see that when the domain Ω ⊂ R d is replaced by λΩ, and u(x) is replaced by u λ (x) := u(x/λ), then holds.Similarly, one can show the following result for (1.2).
Lemma 2.3.Suppose that λ > 0. Then (2.6) Lemma 2.3 shows that the total variation of the normal (1.2) will go to zero when the domain Ω degenerates to a point as λ → 0. This is to be expected since the total variation (1.1) behaves in the same way.In practice, this will not be an issue since (1.2) will always be combined with appropriate data fidelity terms.By contrast (2.5) proposed in Elsey, Esedoglu, 2009 is invariant w.r.t.scaling and thus, in this particular respect, does not generalize (1.1).

Elements of Shape Calculus.
In this section we briefly recall some elements of shape calculus, as necessary in order to study optimization problems in which the domain Ω ⊂ R 3 is an optimization variable.Here we follow common practice and define transformations of Ω in terms of perturbations of identity.That is, we consider families of perturbed domains Ω ε whose material points are given by (2.7) Here V : D → R 3 is some smooth vector field defined on a hold-all D ⊃ Ω. Suppose that J is a functional depending on the domain.Then we are going to denote by dJ(Ω)[V ] the directional shape derivative (also known as Eulerian derivative) of J in the direction of V , i.e., (2.8) Likewise, we write dj(Γ)[V ] for functionals j depending on the surface Γ of Ω.In particular, for an integral of the type the directional shape derivative is given by Delfour, Zolésio, 2011, Equation (4.17) or Schmidt, 2010, Lemma 3.3.13as (2.10) Here •) is the local shape derivative.Moreover, we simply write g instead of g(0, •) in (2.10).
The symbol Dg in (2.10), which we will need occasionally, stands for the "full" derivative (in all three spatial directions) of a function g defined in a neighborhood of Γ.We recall that we are denoting the derivative in tangential directions of functions defined on Γ by the symbol D Γ .Notice that Dg and D Γ g are related by Dg = D Γ g − (Dg)nn .
2.5.Minimizers of the Total Variation of the Normal.When Ω ⊂ R 2 and Γ is a one-dimensional manifold, (1.2) and thus (2.3) reduce to the total absolute curvature Γ |k| ds, where k is the single curvature.It is easy to see that this integral has a minimal value of 2π, which is attained precisely for the boundaries Γ of convex, smoothly bounded domains Ω ⊂ R 2 ; see also Chen, 2000, Chapter 21.1 or Brook, Bruckstein, Kimmel, 2005.
The situation is different for our setting Ω ⊂ R 3 .
Theorem 2.4.Spheres are stationary points for (1.2) among all surfaces Γ of constant area.
Proof.We consider the minimization of (1.2) or equivalently, (2.3), subject to the constraint that the surface area equals the constant A 0 .The Lagrangian associated with this problem is given by Here λ ∈ R is the Lagrange multiplier to be determined below.On the perturbed domain with surface Γ ε with the perturbation according to (2.7), the Lagrangian reads where k 1,ε , k 2,ε are the principal curvatures on Γ ε .This integral is of type (2.9) and its shape derivative at the unperturbed surface, according to (2.10), is given by When Ω is a sphere of radius r, we are going to show that dL(0, λ)[V ] = 0 holds for all perturbation fields V in normal direction and with λ = −1/( √ 2 r).In this setting, the principal curvatures are k 1 (s) = k 2 (s) ≡ 1/r.Consequently, Our strategy now is to show that (2.12) To this end, we utilize that the local shape derivative g [V ] and the material derivative dg[V ] are related by g Sokołowski, Zolésio, 1992, Equation (2.163).We can assume, without loss of generality, that g is extended constant in normal direction.Moreover, due to (2.11), we need to consider only perturbation fields V whose restriction to Γ point in the direction of the normal.Consequently, 0 = (Dg)n = (Dg)V holds.From the chain rule, we thus obtain (2.13) From here, establishing the equality in (2.12) is tedious but straightforward, using integration by parts on Γ.We postpone the details to Appendix C.
As a consequence of (2.12) we obtain the representation for all perturbation fields V parallel to n.Clearly, the term in brackets vanishes identically when λ = −1/( √ 2 r) holds.This concludes the proof.
(1) A numerical study shows that among all ellipsoids of equal area, the sphere is indeed the unique minimizer of (1.3).
(2) The proof utilizes the isotropic nature of (1.2).It continues to hold when the surface area constraint is replaced by a volume constraint, albeit with the different value λ = − √ 2/r 2 .

Split Bregman Iteration
In this section we propose an ADMM iteration which generalizes the split Bregman algorithm for total variation problems proposed in Goldstein, Osher, 2009.As is well known, ADMM methods introduce a splitting of variables so that minimization over individual variables becomes efficient.In our setting, the main variable is the unknown domain Ω.Notice that Ω also determines its boundary Γ as well as its normal vector field n, and we will always consider Γ and n as a function of Ω.
The splitting is achieved through the introduction of a new variable d, which is independent of Ω, Γ and n.At the solution, we require the coupling condition d = D Γ n to hold across Γ.We recall that D Γ n denotes the derivative (pushforward) of n.
Written in terms of Ω and d Notice that for convenience of notation, we represent D Γ n in terms of its actions on the two basis vectors ξ i .
Note also that, at the point s ∈ Γ, the equality d i = (D Γ n) ξ i is one in the tangent space T n(s) S 2 .We therefore introduce Lagrange multipliers λ = (λ 1 , λ 2 ), belonging to the same space, and define the augmented Lagrangian associated with (1.3) as follows, After the usual re-scaling b i := λ i /γ, we can rewrite (3.2) as The main difference to an ADMM method in Euclidean or Hilbert spaces is that the vector fields d i and b i have values in the tangent space T n(•) S 2 .Hence they must be updated whenever Ω and thus the normal vector field n are changing.We now address the individual steps in our ADMM method for (3.3) in detail.

The Shape Optimization
Step.We first address the minimization of (3.3) w.r.t. the shape Ω.Following the standard approach of perturbation of identity, we consider perturbed domains Ω ε as in (2.7) where V : D → R 3 is some smooth vector field defined on a hold-all D ⊃ Ω.We are going to denote by d • [V ] the directional derivatives of scalar quantities in the direction of V , or the material derivatives of functions defined on Γ, respectively.
The derivative of the first term, d (u(Ω), Ω)[V ], is not specified here since it depends on the particular PDE underlying the solution operator u(Ω) and the loss function considered.This derivative can be obtained by standard shape calculus techniques and an example will be given in Section 6.
Next we consider the second term in (3.3).Due to the chosen splitting, the vector fields d i are merely transported along under the perturbation (2.7) and thus we define their perturbed counterparts as As a consequence, their material derivatives vanish and the directional derivative of the term under consideration becomes Here Finally we address the terms The vector fields d i are transported according to (3.4) and thus we need not consider their material derivatives.However, we do need to track the dependencies of (D Γ n) ξ i .Using the product rule yields (3.5) In the last equality we utilized that for any differentiable field F , dF Concerning the first term in (3.5), we notice that dn holds for the material derivatives of the normal vector field n.This result can be found, for instance, in Sokołowski, Zolésio, 1992, Equation (3.168) or Schmidt, 2010, Lemma 3.3.6.Notice that dn[V ] is tangential since Concerning the last term in (3.5), we can show that the material derivatives of ξ i satisfy (3.8) The asymmetry of the formulas in (3.8) stems from the fact that we chose to transport the first basis vector ξ 1 along with (2.7) and then to orthogonalize the second basis vector ξ 2 w.r.t. the first.The details of this derivation are given in Appendix B.
Altogether, we may write the directional derivative of the Augmented Lagrangin This directional derivative is the basis of any shape optimization procedure.After all, the minimization of (3.3) w.r.t. the domain Ω represents a fairly standard shape optimization problem, except for the terms in (3.5).They can be handled, however, conveniently by algorithmic differentiation techniques on the discrete level.In our implementation, we subsequently convert the derivative into a shape gradient by means of an appropriate inner product.The details are given in Section 6.Instead of minimizing (3.3) w.r.t.Ω to a certain accuracy, in practice we only perform one gradient step per iteration.This is in line with Goldstein, Osher, 2009, where a Gauss-Seidel sweep is proposed.

The Total Variation Minimization
Step.Before addressing the minimization of (3.3) w.r.t.d we must note that the data b i at any point s ∈ Γ has to belong to the tangent space T n(s) S 2 .Since the surface Γ and hence the field of normal vectors is changing during the shape optimization step, we must first parallely transport the data b i into the new tangent space.This is achieved via parallel transport along geodesics on S 2 .Suppose for brevity of notation that n − denotes the old normal vector field along the boundary where d 1 , d 2 are sought pointwise in the tangent spaces T n(•) S 2 .We recall that the latter are two-dimensional subspaces of R 3 .At this point it is important to note that the data (D Γ n) ξ i + b i belongs to the same tangent spaces.Therefore, just like in the Euclidean setting, the minimizer of (3.10) can be evaluated explicitly and inexpensively via a pointwise, vectorial shrinkage operation, i.e., Here we abbreviated .

The Multiplier
Notice again that all quantities belong to the subspace T n(•) S 2 of R 3 .
To summarize we outline the split Bregman method for (1.3) as Algorithm 3.1.

5:
Parallely transport the multiplier estimate b along the geodesic from n (k) to n (k+1)

6:
Parallely transport the basis vectors Update the Lagrange multipliers, i.e., set b  The split Bregman method proposed in Section 3 needs to be discretized before it can be numerically realized.More precisely, we need to represent the unknown domain Ω in some discrete fashion.Here we follow a standard approach and represent Ω in terms of a geometrically conforming, simplicial mesh, i.e., consisting of tetrahedra.Consequently, the surface Γ will be given in terms of a geometrically conforming, triangular mesh in R 3 without boundary; see Figure 4.1.
Since the surface Γ is no longer smooth, we cannot directly apply our definition (1.2) of the total variation of the normal to it.Indeed, the normal vector field n is now piecewise constant, and its variation is concentrated in spontaneous changes across edges between triangles, rather than gradual changes expressed by the derivative D Γ n.
To motivate an extension of (1.2) to the discrete setting, we consider a family of smooth approximations Γ ε of a triangulated surface Γ.Some notation is required.In the following, we denote by E an arbitrary edge in the surface mesh Γ.Its Euclidean length is denoted by |E|.Each edge has an arbitrary but fixed orientation, so that its two neighboring triangles can be addressed as T + E and T − E .The corresponding normal vectors, constant on each triangle, are denoted by n + E and denotes the geodesic distance on S 2 , i.e., the angle between the two unit vectors n + E and n − E ; see also Figure 6.2.The family of smooth approximations Γ ε is supposed to be of class C 2 such that the flat triangles are preserved up to a collar of order ε and smoothing occurs in bands of width 2ε around the edges.Such an approximation can be constructed, for instance, by a level-set representation of Γ by means of a signed distance function Φ.Then a family of smooth approximations Γ ε can be obtained as zero level sets of mollifications Φ ϕ ε for sufficiently small ε.Here ϕ ε is the standard Friedrichs mollifier in 3D and denotes convolution.An alternative to this procedure is the Proof.Let us denote the vertices in Γ by V and its edges by E. Since mollification is local, the normal vector is constant in the interior of each triangle minus its collar, which is of order ε.Consequently, changes in the normal vector are confined to a neighborhood of the skeleton.We decompose this area into the disjoint union Here I E,ε are the transition regions around edge E where the normal vector is modified due to mollification, and B V,ε are the regions around vertex V .On I E,ε , we can arrange the basis ξ 1,2 to be aligned and orthogonal to E so that holds, which can be easily evaluated as an iterated integral.In each stripe in I E,ε perpendicular to E, n ε changes monotonically between n + E and n − E , so that the integral along this stripe yields the constant d(n + E , n − E ).Since the length of I E,ε parallel to E is |E| up to terms of order ε, we obtain The contributions to |n ε | T V (Γε) from integration over B V,ε are of order ε since This yields the claim.
Consequently, we will use as an extension of the total variation of the normal functional (1.2) in case of discrete surfaces Γ.We refer to it as the discrete total variation of the normal.

Comparison with Preliminary
Work for Discrete Surfaces.The functional (4.3) has been used previously in the literature.First of all, we mention that it fits into the framework of total variation of manifold-valued functions defined in Giaquinta, Mucci, 2007;Lellmann et al., 2013.Specifically in the context of discrete surfaces, we mention Sullivan, 2005 where H E := |E| Θ E appears as the total mean curvature of the edge E and Θ E is the exterior dihedral angle, which agrees with d(n + E , n − E ), see (4.1).Consequently, (4.3) can be written as E H E .Moreover, (4.3) appears as a regularizer in Wu, Zheng, et al., 2015 within a variational model for mesh denoising but the geodesic distances are approximated for the purpose of numerical solution.
In addition, we are aware of Zhang et al., 2015;Zhong et al., 2018, where was proposed in the context of variational mesh denoising.Notice that in contrast to (4.3), (4.4) utilizes the Euclidean as opposed to the geodesic distance between neighboring normals and is therefore an underestimator for (4.3).
Once again, we are not aware of any work in which (4.3) or its continuous counterpart (1.2) were used as a prior in shape optimization or geometric inverse problems.
4.2.Minimizers of the Discrete Total Variation of the Normal.In this section we investigate some properties of the discrete total variation of the normal.
As can be seen directly from (4.3), Lemma 2.3 continues to hold in the discrete setting, i.e., for all λ > 0. Consequently, when studying minimizers of (4.3), we need to impose a constraint which avoids that Γ degenerates to a point.As in Theorem 2.4, we choose a constraint on the surface area for this purpose and consider the following problem.Given a triangulated surface mesh consisting of vertices V , edges E and facets F , find the mesh with the same connectivity which minimizes To the best of our knowledge, a precise characterization of the minimizers of (4. 5) is an open problem.We conjecture that they depend on the connectivity of the mesh.That is, different triangulations of the same (initial) mesh, e.g., a cube, may yield different minimizers.See also Alexa, Wardetzky, 2011 for a related observation in discrete mean curvature flow.
We do have, however, the following partial result.For the proof, we exploit that (4.3) coincides with the discrete total mean curvature and utilize results from discrete differential geometry.The reader may with to consult Meyer et al., 2003;Polthier, 2005;Wardetzky, 2006;Bobenko, Springborn, 2007;Crane et al., 2013.Theorem 4.2.The icosahedron and the cube with crossed diagonals are stationary for (4.5) among all discrete surfaces Γ of constant area.
Proof.Let us consider the Lagrangian associated with (4.5), Here x i ∈ R 3 denote the coordinates of vertex #i and N V is the total number of vertices of the surface mesh.Notice that the normal vector n ± E , edge lengths |E| and facet areas |F | depend on these coordinates.The gradient of (4.6) w.r.t.x i can be represented as (4.7) see for instance Crane et al., 2013.Here N (i) denotes the index set of vertices adjacent to vertex #i.For any j ∈ N (i), E ij denotes the edge between vertices #i and #j.Moreover, α ij and β ij are the angles as illustrated in Figure 4.3.
For the icosahedron with surface area A 0 , all edges have lengths Moreover, since all facets are unilaterial triangles, α ij = β ij = π/3 holds.Finally, the exterior dihedral angles d(n + Eij , n − Eij ) are all equal to arccos( √ 5/3) ≈ 41.81 • .Consequently, the Lagrangian is stationary for the Lagrange multiplier We remark that (4.3) and thus (4.7) is not differentiable when one or more of the angles d(n + Eij , n − Eij ) are zero.This is the case for the cube with crossed diagonals, see Figure 4.3.However, the right hand side in (4.7) still provides a generalized derivative of L in the sense of Clarke.In contrast to the icosahedron, the cube has two types of vertices.When x i is the center vertex of one of the lateral surfaces, then x i , independently of the value of the Lagrange multiplier λ.Now when x i is a vertex of "corner type", we need to distinguish two types of edges.Along the three edges leading to neighbors of the same type, we have a exterior dihedral angle of Along the three remaining edges leading to surface centers, we have d(n + Eij , n − Eij ) = π/2 and α ij = β ij = π/4.Thus for vertices of "corner type", it is straightforward to verify that 0 belongs to the generalized (partial) differential of L at (x 1 , . . ., x N V , λ) w.r.t.
holds, which is true for the obvious choice of λ.
Numerical experiments show that the icosahedron as well as the cube are not only stationary points, but also local minimizers of (4.5).We can thus conclude that the discrete objective (4.3) exhibits different minimizers than its continuous counterpart (1.2).In particular, (4.3) admits piecewise flat minimizers such as the cube.This property sets our functional apart from other functionals previously used as priors in shape optimization and geometric inverse problems.For instance, the popular surface area functional is well known to produce smooth shapes.
We close this section by comparing the values of (1.2) and (4.3) for the unit cube as well as the for a single regular tetrahedron, for the icosahedron and for (discretized) spheres of the same surface area as the cube.We created triangular meshes of this sphere with various resolutions using Gmsh and evaluated (4.3) numerically.
The results are shown in Tables 4.1 and 4.2.They reveal a factor of approximately √ 2 between the discrete and continuous functionals for the sphere.To explain this discrepancy, notice that the principal curvatures of the sphere are   4.3) for the cube with edge length 1, as well as the regular tetrahedron, the icosahedron and the (discretized) sphere with the same surface area as the cube.
This implies that the derivative map D Γ n has rank two everywhere.Discretized surfaces behave fundamentally different in the following respect.Their curvature is concentrated on the edges, and one of the principal curvatures (in the direction along the edge) is always zero.So even for successively refined meshes, e.g., of the sphere, one is still measuring only one principal curvature at a time.We are thus led to the conjecture that the limit of (4.3) for sucessively refined meshes is the "anisotropic", yet still intrinsic measure Γ |k 1 | + |k 2 | ds, whose value for the sphere in Table 4.1 is 4 √ 6π ≈ 17.3664 and which will be investigated elsewhere.The factor √ 2 can thus be attributed to the ratio between the 1 -and 2 -norms of the vector (1, 1) .Also, one could consider an "isotropic" version of (4.3) in which the dihedral angles across all edges meeting at any given vertex are measured jointly.These alternatives will be considered elsewhere.

Discrete Split Bregman Iteration
In this section we address the discrete realization of the split Bregman iteration presented in Section 3. We continue to write Ω for the unknown domain, but understand that it stands for a tetrahedral mesh whose connectivity is fixed, but whose vertex coordinates will be altered throughout the optimization.The boundary mesh of Ω is a triangulated surface representing the boundary Γ of Ω.We continue to denote the edges of the boundary mesh Γ by E. π , their values of (4.3) and the ration between (4.3) and (1.2).
The discrete analogue of problem (1.3) then reads w.r.t. the vertex positions of Ω. (5.1) As before, u(Ω) denotes the solution of some (discretized) partial differential equation, which depends on the unknown domain Ω.Moreover, represents a loss function.We will consider a concrete example in Section 7.
Notice that the second term in the objective in (5.1) is non-differentiable whenever n + E = n − E occurs on at least one edge.Therefore, similar to Section 3, we introduce a splitting in which the variation of the normal vector becomes an independent variable.Since this variation is confined to edges, where the normal vector jumps (without loss of generality) from n + E to n − E , this new variable becomes Here log n +

E
n − E denotes the logarithmic map, which specifies the unique tangent vector at the point n + E such that the geodesic departing from n + E in that direction will reach n − E at unit time.The logarithmic map is well-defined whenever Together with the set of Lagrange multipliers b E ∈ T n + E S 2 , we define the Augmented Lagrangian pertaining to (5.1), similarly as in (3.3): The vectors d and b are simply the collections of their entries The split Bregman iteration in the discrete setting is very similar to the continuous setting, see Algorithm 3.1.In particular, the TV minimization step can be solved explicitly by one vectorial shrinkage operation per edge E. Specifically, the minimizer of (5.2) with given data Ω (k+1) and associated normal field n (k+1) , as well as multiplier data b 3) for each edge E.
For the sake of completeness, we state the split Bregman iteration for the discrete setting in Algorithm 5.1.

5:
Parallely transport the multiplier estimate b Update the Lagrange multipliers, i.e., set b for all edges E 8: Set k := k + 1 9: end while

Implementation of a Model Problem in FEniCS
In this section we address some details concerning the implementation of Algorithm 5.1 in the finite element framework FEniCS (version 2018.2.dev0), Logg, Mardal, Wells, et al., 2012;Alnaes, Blechta, et al., 2015.For concreteness, we elaborate on a particular loss function (u(Ω), Ω) arising from geological electrical impedance tomography problems with Robin-type far-field boundary conditions.We introduce the problem under consideration first and discuss implementation details and derivative computations later on.
6.1.EIT Model Problem.Electrical impedance tomography (EIT) problems are a prototypical class of inverse problems.Common to these problems is the task of reconstructing the internal conductivity inside a volume from boundary measurements of electric potentials or currents.These problems are both nonlinear and severely ill-posed and require appropriate regularization; see for instance Santosa, Vogelius, 1990;Cheney, Isaacson, Newell, 1999;Chung, Chan, Tai, 2005.Traditionally, EIT problems are modeled with Neumann (current) boundary conditions and the internal conductivity is an unknown function across the entire domain.In order to focus on the demonstration of the utility of (1.2) as a regularizer in geometric inverse problems, we consider a simplified situation in which we seek to reconstruct a perfect conductor inside a domain of otherwise homogeneous electrical properties.
Consequently, the unknown is the interface of the inclusion.As a perfect conductor shields its interior from the electric field, there is no necessity to mesh and simulate the interior of the inclusion.However, we mention that our methodology can be extended also to non-perfect conductors and other geometric inverse problems.The geometry of our model is shown in Figure 6.1, where Γ 1 is the unknown boundary of the perfect conductor and Γ 2 is a fixed boundary where currents are injected and measurements are taken.We assume that i = 1, . . ., r ∈ N experiments are conducted, each resulting in a measured electric potential z i ∈ L 2 (Γ 2 ) on the outer boundary Γ 2 .Experiment #i is conducted by applying the right hand side source f i ∈ L 2 (Γ), which is the characteristic function of one of the colored regions shown in Figure 6.1.We then seek to reconstruct the interface of the inclusion Γ 1 by solving the following regularized least-squares problem of type (1.3), Here u i ∈ H 1 (Ω) is the computed electric field for source f i .Hence, the problem features r PDE constraints with identical operator but different right hand sides.
As detailed in Section 6.2, we compute the derivative of the least-squares objective and the PDE constraint separately from the derivative of the regularization term.
We now focus on the derivative d (u(Ω), Ω) of the data fit term as required for (3.9) in the shape optimization step inside the split Bregman iteration, Algorithm 3.1.The approach is described in the continuous setting but applies similarly to the discrete problem.To evaluate this derivative, we utilize a classical adjoint approach.
To this end, we consider the Lagrangian F (u 1 , . . ., u r , p 1 , . . ., p r Ω) := The differentiation w.r.t.u i leads to the following adjoint problem for p i : (6.3) All states u i and adjoint states p i are discretized in the finite element space CG 1 (Ω) consisting of piecewise linear, globally continuous functions on a tetrahedral mesh covering Ω.Since all forward and adjoint problems are governed by the same differential operator, we assemble the associated stiffness matrix once and solve the state and adjoint equations via an ILU-preconditioned conjugate gradient method.
Provided that u i and p i solve the respective state and adjoint equations, the directional derivative of (u(Ω), Ω) coincides with the partial directional derivative In practice, we evaluate the latter using the coordinate derivative functionality of FEniCS as described in the following subsection.
6.2.Discrete Shape Derivative.Our split Bregman scheme requires the shape derivative of (5.2), which is given by for the problem at hand, where originates from the splitting approach.The term d (u(Ω), Ω)[V ] is computed via the adjoint approach as explained above, Both terms in (6.4) are evaluated using the coordinate derivative capability Ham et al., 2018 in the latest FEniCS release 2018.2.dev0.Being a variation of automatic differentiation (AD), the coordinate derivative provides directional derivatives w.r.t. the vertex position of the mesh over which the respective expression is defined.In case of F , this amounts to the availability of In order to employ this AD functionality, (6.5) needs to be given as a UFL form, a domain specific language based on Python, which forms the native language of the FEniCS framework, see Alnaes, Logg, et al., 2014.Such a UFL representation is easy to achieve if all mathematical expressions are finite element functions.Notice that d and b in (6.5) are constant functions on the edges of the boundary mesh representing Γ 1 .We can thus represent them in the so called HDivTrace space of lowest order in FEniCS.
From the directional derivatives (6.4), we pass to a shape gradient on the surface w.r.t. a scaled H 1 (Γ 1 ) scalar product by solving a variational problem.This problem involves the weak form of a Laplace-Beltrami operator with potential term and it finds W Γ1 ∈ CG 1 (Γ 1 ) 3 such that holds for all test functions V Γ1 ∈ CG 1 (Γ 1 ) 3 .Here P Ω (V Γ1 ) is the extension of V Γ1 to the volume Ω by padding with zeros.
The previous procedure provides us with a shape gradient W Γ1 on the surface Γ 1 alone.In order to propagate this information into the volume Ω, we solve the following mesh deformation equation: find for all test functions V Ω ∈ CG 1 (Ω) 3 with zero Dirichlet boundary conditions, where W Ω is subject to the Dirichlet boundary condition W Ω = W Γ1 and W Ω = 0 on Γ 2 .Subsequently, the vertices of the mesh are moved in the direction of W Ω .
6.3.Intrinsic Formulation Using Co-Normal Vectors.We recall that our functional of interest (4.3) is formulated in terms of the unit outer normal n of the oriented surface Γ 1 .This leads to the term (6.5) inside the augmented Lagrangian (5.2).In order to utilize the differentiation capability of FEniCS w.r.t.vertex coordinates, we need to represent (6.5) in terms of an integral.Since the edges are the interior facets of the surface mesh for Γ 1 , and d and b can be represented as constant on edges as explained above, (6.5) can indeed be written as an integral w.r.t. the interior facet measure dS on Γ 1 .Then, however, the outer normal vectors appearing in the term log n + E n − E is not available.We remedy the situation by observing that the geodesic distance between two normal vectors n + E and n − E on the two triangles T 1 and T 2 sharing the edge E can also be expressed via the conormal (or in-plane normal) vectors µ + E , µ − E , as is shown in Figure 6.2.Indeed, one has Since the co-normal vectors are intrinsic to the surface Γ 1 , they are available on Γ 1 while n + E and n − E are not.

Numerical Results
In this section we present numerical results obtained with Algorithm 5.1 for the geological impedance tomography model problem described in the previous section.The data of the problem are given in Table 7.1.The state u and adjoint state p were discretized using piecewise linear, globally continuous finite elements on a tetrahedral grid of Ω minus the volume enclosed by Γ 1 .The mesh has 4429 vertices and 41 272 tetrahedra.
We show in Figure 7.1 the results obtained in the noise-free setting and with noise.
In the latter case, we added normally distributed random noise with zero mean and standard deviation σ = 10 −2 per degree of freedom on Γ 2 for each of the r = 48 simulations of the forward model (6.1).The amount of noise has to be interpreted in relation to the range of values for the simulated state, which is In each case, the initial guess for Γ 1 was the surface of the ball B 0.5 (0) while the true solution is cube.
T  The results obtained by Algorithm 5.1 applied to the regularized problem (6.1) are shown in Figure 7.1.For comparison, we provide in Figure 7.2 results obtained for a related problem with the same data but with the popular perimeter or surface area regularization, where β |n| T V (Γ1) is replaced by β Γ1 ds = β F |F |.Since in this case the problem is smooth, we applied a shape gradient scheme directly rather than a split Bregman scheme.The regularization parameter β was selected by hand in each case.Automatic parameter selection strategies can clearly be applied here as well but this is out of the scope of the present paper.
As is expected and well known, the use of perimeter regularization leads to results in which the identified inclusion Γ 1 has been smoothed out.This can be explained by the observation that the gradient based minimization of the surface area yields a mean curvature flow.By contrast, our novel prior (4.3) allows for piecewise flat shapes and thus the interface Γ 1 is closely reconstructed in the noise-free situation.Even in the presence of noise the reconstruction is remarkably good.In particular, the flat lateral surfaces and sharp edges can be identified quite well.

Conclusions
In this paper we introduced an analogue of the total-variation prior for the normal vector field.In the continuous setting, this functional is also known as the total root mean square curvature and it admits spheres as minimizers under an area constraint.We also considered a discrete version, which is known as the total mean curvature.While we are currently unable to characterize all of its minimizers, we showed that the icosahedron and a cube with crossed diagonals are stationary under an area constraint.We conjecture that the full set of minimizers is much richer than this.see for instance Hosseini, Uschmajew, 2017 andPersch, 2018, Section 2.3.1, repectively.Here we used the abbreviations v = log n n , |v| g = d(n, n ) and u = v |v|g .To see that both expressions in (A.5) coincide -after plugging in the definition of the geodesic distance (A.1)-it remains to show that which holds true since the norm of the logarithmic map is Hence multiplying with the denominator of the first term in (A.5) yields the equality with the second term, since using the definition of the logarithmic map we obtain

B Material Derivative of the Tangent Basis
In this section we derive (3.8) by a constructive approach.Beginning from a parametrization h of the surface, we give an explicit formula for the tangent basis.The perturbed surface Γ ε is then expressed via a perturbed parameterization h ε := T ε • h, where T ε is given by (2.7).We derive a formula for the perturbed tangent basis via the Gram-Schmidt process.The desired material derivatives are then given by the total derivative w.r.t.ε = 0.
Let h : R 2 → R 3 be an orthogonal parametrization of Γ, i.e. the derivative Dh is a matrix with orthonormal columns, such that s ∈ Γ is locally given by s = h(x) for some x ∈ R 2 .Hence, we define the orthonormal tangent vectors ξ 1 , ξ 2 via Regarding ξ 2 , we do it in a similar way, but have to apply the Gram-Schmidt process to obtain orthonormal perturbed tangent vectors.Hence, ξ 2,ε is given by A straightforward differentiation with respect to ε = 0 results in the material derivatives given in (3.8).
C Details in the Proof of Theorem 2.4 In order to complete the proof of Theorem 2.4, we need to show that Γ 1 g(s) r 2 .To this end, we need a tangential Stokes formula as given in the following lemma.
Lemma C.1.Suppose that a, b are C 1 -vector fields on Γ with values in R 3 , and that V is a C 1 -vector field which is normal, i.e., V = (V n) n holds on Γ.Then we have for all C 1 -vector fields V .We split V in normal and tangential components according to V = (V n) n + 2 i=1 (V ξ i ) ξ i , we arrive at In the last step we used that V is normal and thus V ξ i = 0 holds.
We shall also utilize that g(s) = √ 2/r is a constant on the sphere of radius r.

Figure 2 . 1 .
Figure 2.1.The figure shows part of a smooth surface Γ (blue) and a representation of its tangent spaces at three points s (light blue).The normal vectors are shown in orange.The figure also illustrates that T n(s) S 2 is parallel to T s Γ.
Update.An update of the Lagrange multipliers (b 1 , b 2 ) is achieved, parallel to the Euclidean setting, by replacing b i by b

Figure 4 . 1 .
Figure 4.1.Volume mesh of a cube domain Ω consisting of tetrahedra (left) and corresponding triangular mesh of the boundary Γ (right).

9 :
Set k := k + 1 10: end while 4 Discrete Total Variation of the Normal

Figure 4 . 2 .
Figure 4.2.Illustration of the approximation of a portion of a triangulated surface Γ (left) by a family of smooth surfaces Γ ε (right).Two vertex caps B V,ε and one transition region along an edge I E,ε are highlighted, see the proof of Theorem 4.1.

Figure 4 . 3 .
Figure 4.3.The icosahedron and the cube with crossed diagonals, two stationary surfaces for (4.5).The highlighted regions as well as the figure on the right illustrate the proof of Theorem 4.2.

Figure 6
Figure 6.1.The left plot shows the domain Ω considered in the numerical example.Each color on the outer boundary represents the support of one out of r = 48 electric sources f i .The right figure shows a wireframe plot revealing the true inclusion Γ 1 , i.e., the boundary of the cube.
Figure 6.2.The geodesic distance between normals n + E and n − E (shown in black) of two triangles T + E , T − E which share the edge E agrees with the geodesic distance between the co-normals µ + E and −µ − E (shown in orange).

Figure 7
Figure7.1.Numerical reconstruction of the true (cube shaped) inclusion using Algorithm 5.1 for the solution of (6.1) after 500 iterations in the noise-free situation (left) and after 463 iterations with noise added (right).For the data see Table7.1.

Figure 7 . 2 .
Figure7.2.Numerical reconstruction of the true (cube shaped) inclusion using a shape gradient descent scheme for the solution of (6.1) with perimeter regularization, i.e., with β |n| T V (Γ1) replaced by β Γ1 ds = β F |F |.For the data see Table7.1.The results shown have been obtained after 600 iterations in the noise-free situation (left) and after 600 iterations with noise added (right).
Finally, we utilize (Dn)(s) ≡ id r and (D Γ n)(s) = id r id −nn (C.4) and thus (D Γ n) ξ = ξ/r holds for i = 1, 2. The three terms contributing to the material derivative d[(D Γ n) ξ i ][V ] in (C.1) are given in (3.5) and we consider them individually.We utilize that the Riemannian metric on S 2 is the Euclidean inner product of the ambient R 3 , i.e., g(a, b) = a b.First Term.The insertion of the first term in (3.5) into the left hand side of (CΓ (dn[V ]) ξ i ds by (C.4) = 1 √ 2 Γ div Γ dn[V ] ds = 0. (C.5) the geodesic from n − to n.This step is inexpensive since geodesics and parallel transport on S 2 are available in terms of explicit formulas; see Appendix A.We can now address the minimization of (3.3) w.r.t.d = (d 1 , d 2 ).Since the first term in (3.3) does not depend on d, we are left with the minimization of

Table 4
.1.Comparison between the continuous and the discrete total variation of the normal functionals (1.2) and (

Table 7 .
1. Setting of the numerical experiments.