Total Variation of the Normal Vector Field as Shape Prior

TL;DR

This paper introduces a total variation prior for the normal vector field on 3D shape boundaries, analyzing its properties using differential geometry and proposing an extension of the split Bregman method for non-differentiable cases.

Contribution

It develops a novel total variation prior for normal vector fields on 3D shapes and extends optimization methods to handle non-differentiability on flat regions.

Findings

Spheres are stationary points under area constraints.

The total variation functional is non-differentiable on flat boundary regions.

An extended split Bregman method for manifold-valued functions is proposed.

Abstract

An analogue of the total variation prior for the normal vector field along the boundary of smooth shapes in 3D is introduced. The analysis of the total variation of the normal vector field is based on a differential geometric setting in which the unit normal vector is viewed as an element of the two-dimensional sphere manifold. It is shown that spheres are stationary points when the total variation of the normal is minimized under an area constraint. Shape calculus is used to characterize the relevant derivatives. Since the total variation functional is non-differentiable whenever the boundary contains flat regions, an extension of the split Bregman method to manifold valued functions is proposed.

Equations

$$\lvert u\rvert_{TV(\Omega)}\coloneqq\int_{\Omega}\lvert\nabla u\rvert_{2}\,\textup{d}x=\int_{\Omega}\mathopen{\bigl{(}}{\lvert(Du)\,{\boldsymbol{e}}_{1}\rvert^{2}+\lvert(Du)\,{\boldsymbol{e}}_{2}\rvert^{2}}\mathclose{\bigr{)}}.$$

$$\lvert{\boldsymbol{n}}\rvert_{TV(\Gamma)}\coloneqq\int_{\Gamma}\mathopen{\bigl{(}}{\lvert(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{1}\rvert_{\mathfrak{g}}^{2}+\lvert(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{2}\rvert_{\mathfrak{g}}^{2}}\mathclose{\bigr{)}}^{1/2}\,\textup{d}{\boldsymbol{s}}$$

$$\int_{\Gamma}\mathopen{\bigl{(}}{\lvert(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{1}\rvert_{\mathfrak{g}}^{2}+\lvert(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{2}\rvert_{\mathfrak{g}}^{2}}\mathclose{\bigr{)}}^{1/2}\,\textup{d}{\boldsymbol{s}}=\int_{\Omega}\mathopen{\bigl{(}}{\textstyle\bigl{\lvert}\frac{\partial u}{\partial x_{1}}\bigr{\rvert}^{2}+\bigl{\lvert}\frac{\partial u}{\partial x_{2}}\bigr{\rvert}^{2}}\mathclose{\bigr{)}}^{1/2}\,\textup{d}x=\displaystyle\int_{\Omega}\lvert\nabla u\rvert_{2}\,\textup{d}x.$$

$$\text{Minimize}\ell (u(\Omega ),\Omega )+\beta |\mathbf{n}{|}_{TV(\Gamma )}$$

$$\text{w.r.t. Ω in a suitable class of domains}.$$

$$\begin{gathered} |(D_{\Gamma }\mathbf{n}){\mathbf{\xi }}_1{|}_{\mathfrak{g}}^2+|(D_{\Gamma }\mathbf{n}){\mathbf{\xi }}_2{|}_{\mathfrak{g}}^2=\mathrm{trace}(J^{\top }J)=\mathrm{trace}(J^{\top }JQ{Q}^{\top })=\mathrm{trace}(Q^{\top }{J}^{\top }JQ) \\ =|(D_{\Gamma }\mathbf{n})Q{\mathbf{\xi }}_1{|}_{\mathfrak{g}}^2+|(D_{\Gamma }\mathbf{n})Q{\mathbf{\xi }}_2{|}_{\mathfrak{g}}^2=|(D_{\Gamma }\mathbf{n}){\mathbf{\eta }}_1{|}_{\mathfrak{g}}^2+|(D_{\Gamma }\mathbf{n}){\mathbf{\eta }}_2{|}_{\mathfrak{g}}^2. \end{gathered}$$

$$S:{\mathcal{T}}_{\mathbf{s}}\Gamma \to {\mathcal{T}}_{\mathbf{s}}\Gamma .$$

$$\mathopen{\bigl{(}}{\lvert(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{1}\rvert_{\mathfrak{g}}^{2}+\lvert(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{2}\rvert_{\mathfrak{g}}^{2}}\mathclose{\bigr{)}}^{1/2}=\mathopen{\bigl{(}}{k_{1}^{2}+k_{2}^{2}}\mathclose{\bigr{)}}^{1/2}.$$

$$|(D_{\Gamma }\mathbf{n}){\mathbf{\xi }}_1{|}_{\mathfrak{g}}^2+|(D_{\Gamma }\mathbf{n}){\mathbf{\xi }}_2{|}_{\mathfrak{g}}^2=(S{\mathbf{\xi }}_1{)}^{\top }(S{\mathbf{\xi }}_1)+(S{\mathbf{\xi }}_2{)}^{\top }(S{\mathbf{\xi }}_2).$$

$$|(D_{\Gamma }\mathbf{n}){\mathbf{\xi }}_1{|}_{\mathfrak{g}}^2+|(D_{\Gamma }\mathbf{n}){\mathbf{\xi }}_2{|}_{\mathfrak{g}}^2=k_1^2|{\mathbf{\xi }}_1{|}_{\mathfrak{g}}^2+k_2^2|{\mathbf{\xi }}_2{|}_{\mathfrak{g}}^2=k_1^2+k_2^2.$$

$$|\mathbf{n}{|}_{TV(\Omega )}={\int }_{\Gamma }(k_1^2+k_2^2{)}^{1/2}\text{d}\mathbf{s},$$

$${\int }_{\Gamma }(k_1^2+k_2^2)\text{d}\mathbf{s}$$

$${\int }_{\Gamma }|k_1{k}_2|\text{d}\mathbf{s}$$

$$|{\mathbf{n}}_{\delta }{|}_{TV(\delta \Gamma )}=\delta |\mathbf{n}{|}_{TV(\Gamma )}.$$

$${\mathbf{x}}_{\epsilon }={\mathbf{T}}_{\epsilon }(\mathbf{x}):=\mathbf{x}+\epsilon \mathbf{V}(\mathbf{x}).$$

$$\text{d}J(\Omega )[\mathbf{V}]=\underset{\epsilon \searrow 0}{\lim }\frac{J({\Omega }_{\epsilon })-J(\Omega )}{\epsilon }.$$

$$J({\Gamma }_{\epsilon })={\int }_{{\Gamma }_{\epsilon }}g(\epsilon ,{\mathbf{s}}_{\epsilon })\text{d}{\mathbf{s}}_{\epsilon },$$

$$\text{d}J(\Gamma )[\mathbf{V}]={\int }_{\Gamma }g(0,\mathbf{s}){\mathrm{div}}_{\Gamma }\mathbf{V}(\mathbf{s})+\text{d}g[\mathbf{V}](0,\mathbf{s})\text{d}\mathbf{s},$$

$$\text{d}g[\mathbf{V}](0,\mathbf{s})={\frac{\text{d}}{\text{d}\epsilon }|}_{\epsilon =0}g(\epsilon ,{\mathbf{s}}_{\epsilon })={\frac{\text{d}}{\text{d}\epsilon }|}_{\epsilon =0}g(\epsilon ,{\mathbf{T}}_{\epsilon }(\mathbf{s}))$$

$${\mathrm{div}}_{\Gamma }\mathbf{V}=\sum _{i=1}^2{\mathbf{\xi }}_i^{\top }(D\mathbf{V}){\mathbf{\xi }}_i=\mathrm{div}\mathbf{V}-{\mathbf{n}}^{\top }(D\mathbf{V})\mathbf{n}.$$

$$g^{\prime }[\mathbf{V}]=\text{d}g[\mathbf{V}]-(Dg)\mathbf{V}.$$

$$\textup{d}J(\Gamma)[{\boldsymbol{V}}]=\int_{\Gamma}{\boldsymbol{V}}^{\top}{\boldsymbol{n}}\mathopen{\bigl{[}}{(Dg){\boldsymbol{n}}+(k_{1}+k_{2})\,g}\mathclose{\bigr{]}}+g^{\prime}[{\boldsymbol{V}}]\,\textup{d}{\boldsymbol{s}},$$

$$D(\text{d}\mathbf{F}[\mathbf{V}])$$

$$=(D\mathbf{F}{)}^{\prime }[\mathbf{V}]+D((D\mathbf{F})\mathbf{V})$$

$$=(D\mathbf{F}{)}^{\prime }[\mathbf{V}]+(D(D\mathbf{F})\mathbf{V})+(D\mathbf{F})(D\mathbf{V})$$

$$=\text{d}(D\mathbf{F})[\mathbf{V}]+(D\mathbf{F})(D\mathbf{V}).$$

$$\int_{\Gamma}{\boldsymbol{a}}^{\top}(D_{\Gamma}{\boldsymbol{V}})\,{\boldsymbol{b}}\,\textup{d}{\boldsymbol{s}}\\ =\int_{\Gamma}{\boldsymbol{V}}^{\top}{\boldsymbol{n}}\mathopen{\bigl{[}}{-\operatorname{div}_{\Gamma}(({\boldsymbol{a}}^{\top}{\boldsymbol{n}})\,{\boldsymbol{b}})+({\boldsymbol{a}}^{\top}{\boldsymbol{n}})({\boldsymbol{b}}^{\top}{\boldsymbol{n}})\,(k_{1}+k_{2})+{\boldsymbol{a}}^{\top}(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{b}}}\mathclose{\bigr{]}}\,\textup{d}{\boldsymbol{s}}.$$

$${\int }_{\Gamma }c{\mathrm{div}}_{\Gamma }\mathbf{V}\text{d}\mathbf{s}={\int }_{\Gamma }{\mathbf{V}}^{\top }\mathbf{n}c(k_1+k_2)\text{d}\mathbf{s}-{\int }_{\Gamma }(D_{\Gamma }c)\mathbf{V}\text{d}\mathbf{s}$$

$${\int }_{\Gamma }{\mathbf{a}}^{\top }(D_{\Gamma }\mathbf{V})\mathbf{b}\text{d}\mathbf{s}={\int }_{\Gamma }{\mathbf{a}}^{\top }{D}_{\Gamma }(({\mathbf{V}}^{\top }\mathbf{n})\mathbf{n})\mathbf{b}+\sum _{i=1}^2{\mathbf{a}}^{\top }{D}_{\Gamma }(({\mathbf{V}}^{\top }{\mathbf{\xi }}_i){\mathbf{\xi }}_i)\mathbf{b}\text{d}\mathbf{s}$$

$$={\int }_{\Gamma }{D}_{\Gamma }({\mathbf{V}}^{\top }\mathbf{n})({\mathbf{a}}^{\top }\mathbf{n})\mathbf{b}+({\mathbf{V}}^{\top }\mathbf{n}){\mathbf{a}}^{\top }(D_{\Gamma }\mathbf{n})\mathbf{b}$$

$$+\sum _{i=1}^2{D}_{\Gamma }({\mathbf{V}}^{\top }{\mathbf{\xi }}_i)({\mathbf{a}}^{\top }{\mathbf{\xi }}_i)\mathbf{b}+({\mathbf{V}}^{\top }{\mathbf{\xi }}_i){\mathbf{a}}^{\top }(D_{\Gamma }{\mathbf{\xi }}_i)\mathbf{b}\text{d}\mathbf{s}\text{(by the product rule)}$$

$$={\int }_{\Gamma }{\mathbf{V}}^{\top }\mathbf{n}\left[({\mathbf{a}}^{\top }\mathbf{n})({\mathbf{n}}^{\top }\mathbf{b})(k_1+k_2)-{\mathrm{div}}_{\Gamma }(({\mathbf{a}}^{\top }\mathbf{n})\mathbf{b})+{\mathbf{a}}^{\top }(D_{\Gamma }\mathbf{n})\mathbf{b}\right]$$

$$+\sum _{i=1}^2{\mathbf{V}}^{\top }{\mathbf{\xi }}_i\left[{\mathbf{a}}^{\top }(D_{\Gamma }{\mathbf{\xi }}_i)\mathbf{b}-{\mathrm{div}}_{\Gamma }(({\mathbf{a}}^{\top }{\mathbf{\xi }}_i)\mathbf{b})\right]\text{d}\mathbf{s}(\text{by eqrefeq:TangentStokesVanilla})$$

$$={\int }_{\Gamma }{\mathbf{V}}^{\top }\mathbf{n}\left[({\mathbf{a}}^{\top }\mathbf{n})({\mathbf{n}}^{\top }\mathbf{b})(k_1+k_2)-{\mathrm{div}}_{\Gamma }(({\mathbf{a}}^{\top }\mathbf{n})\mathbf{b})+{\mathbf{a}}^{\top }(D_{\Gamma }\mathbf{n})\mathbf{b}\right].$$

$$\text{d}\mathbf{n}[\mathbf{V}]=-(D_{\Gamma }\mathbf{V}{)}^{\top }\mathbf{n}.$$

Full text

1   Introduction

The total variation (TV) functional is popular as a regularizer in imaging and inverse problems; see for instance [RudinOsherFatemi1992, ChanGolubMulet1999, BachmayrBurger2009, Langer2017] and [Vogel2002, Chapter 8]. For a real-valued function $u\in W^{1,1}(\Omega)$ on a bounded domain $\Omega$ in $\mathbb{R}^{2}$, the total variation seminorm is defined as

[TABLE]

Notice that we restrict the discussion to the isotropic case here, i.e., $\lvert\,\cdot\,\rvert_{2}$ denotes the Euclidean norm. Moreover, $Du$ is the derivative of $u$ and $\{{\boldsymbol{e}}_{1},{\boldsymbol{e}}_{2}\}$ denotes the standard Euclidean basis. The seminorm (1.1) extends to less regular, so-called BV functions (bounded variation), whose distributional gradient exists only in the sense of measures. We refer the reader to [Giusti1984, AttouchButtazzoMichaille2006] for an extensive discussion of BV functions. 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, which can be used, for instance, in shape optimization applications as well as geometric inverse problems. In the latter class, the unknown, which one seeks to recover, is a shape $\Omega\subset\mathbb{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 $\Omega$ will be denoted by $\Gamma$.

The novel functional, which we term the total variation of the normal field along a smooth surface $\Gamma$, is defined by

[TABLE]

in analogy to (1.1). In (1.2), ${\boldsymbol{n}}$ is the outer unit normal vector field along $\Gamma$, i.e., ${\boldsymbol{n}}$ belongs to the manifold $\SS^{2}=\{{\boldsymbol{v}}\in\mathbb{R}^{3}:\lvert{\boldsymbol{v}}\rvert_{2}=1\}$ pointwise. Moreover, $D_{\Gamma}{\boldsymbol{n}}$ denotes the derivative (push-forward) of the normal vector field, and $\{{\boldsymbol{\xi}}_{1}({\boldsymbol{s}}),{\boldsymbol{\xi}}_{2}({\boldsymbol{s}})\}$ denotes an orthonormal basis (w.r.t. the Euclidean inner product in the embedding $\Gamma\subset\mathbb{R}^{3}$) of the tangent spaces $\mathcal{T}_{{\boldsymbol{s}}}\Gamma$ along $\Gamma$. Finally, $\lvert\,\cdot\,\rvert_{\mathfrak{g}}$ denotes the norm induced by a Riemannian metric on $\SS^{2}$. We will consider the metric induced from embedding $\SS^{2}$ in $\mathbb{R}^{3}$, i. e., the distance induced by this metric is the arc length distance and the curvature is 1. We write $\lvert\,\cdot\,\rvert_{\mathfrak{g}}$ for the norm induced by the Riemannian metric $\mathfrak{g}(\cdot,\,\cdot)$ on the tangent spaces $\SS^{2}$.

Let us argue that (1.2) generalizes (1.1). Since the normal vector field ${\boldsymbol{n}}$ replaces the scalar-valued function $u$ in (1.1), assume for the moment that ${\boldsymbol{n}}$ maps into $\mathbb{R}$ instead of $\SS^{2}$. Then the tangent space $\mathcal{T}_{{\boldsymbol{n}}}\mathbb{R}$ is equal to $\mathbb{R}$, endowed with its usual inner product. Finally, the manifold $\Gamma$ in (1.2) takes the role of $\Omega\subset\mathbb{R}^{2}$ in (1.1). We can choose, without loss of generality, the basis ${\boldsymbol{\xi}}_{i}={\boldsymbol{e}}_{i}$. Consequently, (1.2) becomes (1.1):

[TABLE]

A thorough introduction to $\lvert{\boldsymbol{n}}\rvert_{TV(\Gamma)}$ 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. (1)

The variable on which (1.2) depends is the domain $\Omega$. Since the normal vector field ${\boldsymbol{n}}$ in turn depends on $\Omega$, both the integration domain $\Gamma$ and the integrand in (1.2) depend on $\Omega$. By contrast, $\Omega$ is fixed in (1.1), where $u$ is the variable. 2. (2)

The normal vector field, whose pointwise variation the total variation functional (1.2) seeks to capture, is manifold-valued with values in $\SS^{2}$. By contrast, the function $u$ in (1.1) is real-valued. 3. (3)

It is well known that the TV functional penalizes jumps and non-zero gradients of BV functions. Consequently, the minimization of (1.1) avoids unnecessary variations of $u$ and thus favors piecewise constant minimizers in BV. Generally, it does not admit minimizers in spaces of functions of higher smoothness, such as $W^{1,1}(\Omega)$. The situation is slightly different for (1.2) since we are considering closed surfaces $\Gamma$, which yields a periodicity constraint for the normal vector field ${\boldsymbol{n}}$. In this setting, unnecessary variations of ${\boldsymbol{n}}$ correspond to non-convex regions of the enclosed body $\Omega$. Consequently, the minimization of (1.2) favors convex shapes and, more precisely, spheres; see Theorem 3.6.

Further properties of (1.2) will be discussed in Section 2.

In this paper we are also considering the total variation of the normal (1.2) as a prior in shape optimization problems, which may involve a partial differential equation (PDE). The aforementioned problem can be cast in the form

[TABLE]

Here $u(\Omega)$ denotes the solution of the problem specific PDE, which depends on the unknown domain $\Omega$. Moreover, $\ell$ represents a loss function, such as a least squares function.

The coupling between $\Omega$ and its normal vector field ${\boldsymbol{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 $\Gamma$, (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 [GoldsteinOsher2009], where ${\boldsymbol{d}}=\nabla u$ was introduced as an independent variable in the context of the total variation functional (1.1). The variables $u$ and ${\boldsymbol{d}}$ are coupled through a constraint, which is then handled in an Alternating Direction Method of Multipliers (ADMM) framework. We refer the reader to [GlowinskiMarroco1975, GoldsteinBressonOsher2010, GoldsteinODonoghueSetzerBaraniuk2014] for more on ADMM.

In our proposed splitting, we introduce a new variable ${\boldsymbol{d}}$, independent of $\Omega$ and its normal vector field ${\boldsymbol{n}}$, and require the coupling condition ${\boldsymbol{d}}=D_{\Gamma}{\boldsymbol{n}}$ to hold across $\Gamma$. An outstanding feature of the proposed splitting is that the two subproblems, the minimization w.r.t. $\Omega$ and w.r.t. ${\boldsymbol{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.

Although many optimization algorithms have recently been generalized to Riemannian manifolds [Bacak2014, BergmannPerschSteidl2016, BergmannHerzogTenbrinckVidal-Nunez2019_preprint], 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 [AbsilMahonySepulchre2008].

The structure of the paper is as follows. In Section 2 we provide an analysis of (1.2) and its properties. We also compare (1.2) to geometric functionals appearing elsewhere in the literature. In section 3 we discuss the role of (1.2) in optimization problems. Section 4 is devoted to the formulation of an ADMM method which generalizes the split Bregman algorithm to the manifold-valued problem (1.3).

This paper is accompanied by a companion paper [BergmannHerrmannHerzogSchmidtVidalNunez2019:2_preprint] where we introduce a discrete counterpart of (1.2) on simplicial meshes, as they are frequently used in finite element discretizations of PDEs. On these piecewise flat surfaces, the normal vector ${\boldsymbol{n}}$ jumps across edges and the definition (1.2) needs to be generalized. An appropriate discrete version of the split Bregman iteration will be presented in the companion paper, along with numerical results for geometric inverse problems.

2   Total Variation of the Normal

In this section we discuss our proposal (1.2) for the total variation of the normal on smooth surfaces 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 [DoCarmo1976, GrayAbbenaSalamon2006, Kuehnel2013] for a thorough introduction.

2.1.  Preliminaries

From this section onwards we assume that the boundary $\Gamma$ of the unknown bounded domain $\Omega$ is a smooth, compact, orientable manifold of dimension 2 without boundary, embedded in $\mathbb{R}^{3}$. Therefore we can think of tangent vectors at ${\boldsymbol{s}}\in\Gamma$ to be elements of the appropriate two-dimensional subspace (the tangent plane) of $\mathbb{R}^{3}$. This tangent plane at ${\boldsymbol{s}}$ is denoted by $\mathcal{T}_{{\boldsymbol{s}}}\Gamma$. Each tangent plane is endowed with the Riemannian metric furnished by the embedding via the pull-back of the Euclidean metric in $\mathbb{R}^{3}$. In other words, the inner product of two vectors ${\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2}\in\mathcal{T}_{{\boldsymbol{s}}}\Gamma$ is simply given by $\mathfrak{g}({\boldsymbol{\xi}}_{1},\,{\boldsymbol{\xi}}_{2})={\boldsymbol{\xi}}_{1}^{\top}{\boldsymbol{\xi}}_{2}$. In what follows, $\{{\boldsymbol{\xi}}_{1}({\boldsymbol{s}}),{\boldsymbol{\xi}}_{2}({\boldsymbol{s}})\}$ denotes an orthonormal basis in $\mathcal{T}_{{\boldsymbol{s}}}\Gamma$. As the following remark shows, the choice of this basis and how it varies with ${\boldsymbol{s}}$ will not matter.

Outward pointing unit normal vectors ${\boldsymbol{n}}$ along $\Gamma$ will be considered elements of the two-dimensional smooth manifold $\SS^{2}$. The derivative or push-forward of the normal map ${\boldsymbol{n}}$ is denoted by $D_{\Gamma}{\boldsymbol{n}}$. At a given ${\boldsymbol{s}}\in\Gamma$, $D_{\Gamma}{\boldsymbol{n}}$ thus maps tangent vectors ${\boldsymbol{\xi}}\in\mathcal{T}_{{\boldsymbol{s}}}\Gamma$ into tangent vectors $(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}\in\mathcal{T}_{{\boldsymbol{n}}({\boldsymbol{s}})}\SS^{2}$. In what follows, we will suppress the dependence on the point ${\boldsymbol{s}}\in\Gamma$ where possible.

Remark 2.1**.**

The total variation of the normal (1.2) is independent of the choice of the orthonormal basis in the tangent spaces $\mathcal{T}_{{\boldsymbol{s}}}\Gamma$. To show this, it is enough to consider a point ${\boldsymbol{s}}\in\Gamma$ and suppose that $\{{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2}\}$ and $\{{\boldsymbol{\eta}}_{1},{\boldsymbol{\eta}}_{2}\}$ are two orthonormal bases of $\mathcal{T}_{{\boldsymbol{s}}}\Gamma$. Then there exists an orthogonal matrix $Q\in\mathbb{R}^{3\times 3}$ such that ${\boldsymbol{\eta}}_{i}=Q\,{\boldsymbol{\xi}}_{i}$ holds for $i=1,2$. For $J\coloneqq\begin{bmatrix}(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{1}&(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{2}\end{bmatrix}$ the integrand in (1.2) satisfies

[TABLE]

Similarly, as we do for $\Gamma$, we consider $\SS^{2}$ embedded into $\mathbb{R}^{3}$ and therefore we can conceive the tangent space $\mathcal{T}_{{\boldsymbol{n}}({\boldsymbol{s}})}\SS^{2}$ as a two-dimensional plane in $\mathbb{R}^{3}$ tangent to the sphere $\SS^{2}$. We endow $\mathcal{T}_{{\boldsymbol{n}}({\boldsymbol{s}})}\SS^{2}$ with the Riemannian metric furnished by the pull-back of the Euclidean metric as well, which we denote by $\mathfrak{g}(\cdot,\,\cdot)$ to distinguish it from the Riemannian metric on $\mathcal{T}_{{\boldsymbol{s}}}\Gamma$. In fact, $\mathcal{T}_{{\boldsymbol{n}}({\boldsymbol{s}})}\SS^{2}$ is clearly parallel to $\mathcal{T}_{{\boldsymbol{s}}}\Gamma$, see Figure 2.1. We can therefore identify the two tangent spaces and we write $\mathcal{T}_{{\boldsymbol{n}}({\boldsymbol{s}})}\SS^{2}\cong\mathcal{T}_{{\boldsymbol{s}}}\Gamma$ to indicate this.

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 ${{\boldsymbol{N}}_{\Gamma}}:\Gamma\to\SS^{2}$ is also known as the Gauss map; see for instance [Kuehnel2013, Chapter 3]. Its derivative at ${\boldsymbol{s}}\in\Gamma$ maps tangent directions in $\mathcal{T}_{{\boldsymbol{s}}}\Gamma$ into tangent directions in $\mathcal{T}_{{\boldsymbol{n}}({\boldsymbol{s}})}\SS^{2}\cong\mathcal{T}_{{\boldsymbol{s}}}\Gamma$. With the latter identification, the derivative of the Gauss map is known as the shape operator

[TABLE]

Notice that $S$ is self-adjoint, i.e., $(S{\boldsymbol{\xi}}_{1})^{\top}{\boldsymbol{\xi}}_{2}=(S{\boldsymbol{\xi}}_{2})^{\top}{\boldsymbol{\xi}}_{1}$ holds for all ${\boldsymbol{s}}\in\Gamma$ and all ${\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2}\in\mathcal{T}_{{\boldsymbol{s}}}\Gamma$; see for instance [GrayAbbenaSalamon2006, Lemma 13.14]. The two eigenvalues of $S$ are the principal curvatures of the surface $\Gamma$ at ${\boldsymbol{s}}$, denoted by $k_{1}$ and $k_{2}$. This insight allows us to interpret the integrand in (1.2) differently.

Proposition 2.2**.**

The integrand in (1.2) satisfies

[TABLE]

Proof.

Consider the square of the integrand,

[TABLE]

Due to Remark 2.1 we can choose ${\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2}$ to be normalized eigenvectors in $\mathcal{T}_{{\boldsymbol{s}}}\Gamma$ corresponding to the eigenvalues $k_{1}$, $k_{2}$, respectively. Therefore we get

[TABLE]

∎

2.3.  Comparison with Prior Work

The representation of the integrand from Proposition 2.2 allows us to rewrite (1.2) as the integral over the root mean square curvature,

[TABLE]

and compare it with related functionals appearing in the literature. The quantity

[TABLE]

is known as the integral over the total curvature (although this term is also used for other quantities in the literature). The functional (2.3) has a long tradition in surface fairing applications and can be interpreted as a surface strain energy, see for instance [LottPullin1988, HagenSchulze1987, WelchWitkin1992, HalsteadKassDeRose1993, WelchWitkin1994, Greiner1994, TasdizenWhitakerBurchardOsher2003]. Since (2.3) corresponds to $\int_{\Omega}\lvert\nabla u\rvert_{2}^{2}\,\textup{d}x$ in imaging applications, which leads to a Laplacian in the associated optimality conditions (and thus also in the corresponding $L^{2}$-gradient flow), (2.3) tends to smooth the surface and its features. We also mention [KimmelSochen2002] where smoothing by diffusion was employed to $\SS^{1}$-valued images via a mean curvature flow.

By contrast, the functional (2.2) seems to have made very few appearances in the mathematical literature. We are aware of the PhD thesis [Maekawa1993, Chapter 6] and the subsequent book publication [PatrikalakisMaekawa2001] where it was used to guide mesh generation. In [AteshianRosenwasserMow1992, MarzkeTocheriMarzkeFemiani2012] the pointwise root mean square curvature is used as a measure of flatness in biomedical classification problems, in [PullaRazdanFarin2001] for the purpose of surface segmentation and in [WuMaMaJiaHu2010] 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 [McKinnonSpackmanMitchell2004]. We are however not aware of any use of (1.2) or its equivalent form (2.2) as a prior in shape optimization problems.

We regard (1.2) as a natural extension of the total variation seminorm (1.1) to the normal vector on surfaces, measuring surface flatness, but other extensions are certainly possible. Notably, the authors in [ElseyEsedoglu2009] propose the total absolute Gaussian curvature

[TABLE]

for the same purpose. From the Gauss–Bonnet theorem (see for instance [GrayAbbenaSalamon2006, Chapter 27] or [Kuehnel2013, Chapter 4F]) it follows that the boundaries $\Gamma$ of convex domains $\Omega$ will be the global minimizers of (2.4), and they yield a value of $4\pi$. Thus (2.4) 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 $\Omega\subset\mathbb{R}^{d}$ is replaced by ${\delta}\Omega$, and $u(x)$ is replaced by $u_{\delta}(x)\coloneqq u(x/{\delta})$, then $\lvert u_{\delta}\rvert_{TV({\delta}\Omega)}={\delta}^{d-1}\lvert u\rvert_{TV(\Omega)}$ holds. Similarly, we can show that the total variation on a 2-dimensional surface $\Gamma$ scales as follows.

Lemma 2.3**.**

Suppose that ${\delta}>0$. Then

[TABLE]

Lemma 2.3 implies that the total variation of the normal (1.2) will go to zero when the domain $\Omega$ degenerates to a point as ${\delta}\to 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 other, e.g., data fidelity terms. By contrast (2.4) proposed in [ElseyEsedoglu2009] is invariant w.r.t. scaling and thus, in this particular respect, does not generalize (1.1).

3   Analysis of the Total Variation of the Normal

In this section we discuss some properties of the total variation of the normal functional (1.2). To this end, we begin by briefly recalling some elements of shape calculus, as necessary in order to study optimization problems in which the domain $\Omega\subset\mathbb{R}^{3}$ appears as an optimization variable. Then we discuss properties of (1.2). In section 3.3 we briefly comment on the case of curves, i.e., when $\Omega\subset\mathbb{R}^{2}$ and its boundary $\Gamma$ is a one-dimensional manifold.

3.1.  Elements of Shape Calculus

Here we follow common practice and define transformations of $\Omega$ in terms of perturbations of identity. That is, we consider families of perturbed domains $\Omega_{\varepsilon}$ whose material points are given by

[TABLE]

Here ${\boldsymbol{V}}:\mathcal{D}\to\mathbb{R}^{3}$ is some smooth vector field defined on a hold-all $\mathcal{D}\supset\Omega$. Suppose that $J$ is a functional depending on the domain. Then we denote by $\textup{d}J(\Omega)[{\boldsymbol{V}}]$ the directional shape derivative (also known as Eulerian derivative) of $J$ in the direction of ${\boldsymbol{V}}$, i.e.,

[TABLE]

Likewise, we write $\textup{d}J(\Gamma)[{\boldsymbol{V}}]$ for functionals $J$ depending on the surface $\Gamma$ of $\Omega$. In particular, for an integral of the type

[TABLE]

the directional shape derivative is given by [SokolowskiZolesio1992, Eq. (2.172)]

[TABLE]

where the material derivative $\textup{d}g[{\boldsymbol{V}}]$ is defined as the total derivative

[TABLE]

and $\operatorname{div}_{\Gamma}{\boldsymbol{V}}$ denotes the (tangential) divergence of ${\boldsymbol{V}}$ along $\Gamma$. It is related to the divergence in $\mathbb{R}^{3}$ via

[TABLE]

In the following, we simply write $g$ instead of $g(0,\cdot)$ and in addition to the material derivative, we also introduce the (local) shape derivative as the partial derivative $g^{\prime}[{\boldsymbol{V}}]\coloneqq(\partial/\partial\varepsilon)|_{\varepsilon=0}\,g(\varepsilon,\cdot)$. Hence, both are related to each other via

[TABLE]

See for instance [SokolowskiZolesio1992, Eq. (2.163)]. If the shape derivative $g^{\prime}[{\boldsymbol{V}}]$ exists, the shape derivative $\textup{d}J(\Gamma)[{\boldsymbol{V}}]$ can alternatively be expressed as

[TABLE]

see [SokolowskiZolesio1992, Eq. (2.174)]. Furthermore, spatial derivatives and material derivatives of differentiable fields ${\boldsymbol{F}}$ fulfill

[TABLE]

The symbol $Dg$ in (3.5), which we will need occasionally, stands for the “full” derivative (in all three spatial directions) of a function $g$ defined in a neighborhood of $\Gamma$. We recall that we are denoting the derivative in tangential directions of functions defined on $\Gamma$ by the symbol $D_{\Gamma}$. Notice that $Dg$ and $D_{\Gamma}g$ are related by $Dg=D_{\Gamma}g+(Dg){\boldsymbol{n}}{\boldsymbol{n}}^{\top}$.

Lemma 3.1**.**

Suppose that ${\boldsymbol{a}}$, ${\boldsymbol{b}}$ are $C^{1}$-vector fields on $\Gamma$ with values in $\mathbb{R}^{3}$, and that ${\boldsymbol{V}}$ is a $C^{1}$-vector field which is normal, i.e., ${\boldsymbol{V}}=({\boldsymbol{V}}^{\top}{\boldsymbol{n}})\,{\boldsymbol{n}}$ holds on $\Gamma$. Then we have

[TABLE]

Proof.

The general tangential Stokes formula [DelfourZolesio2011, Eq. (5.27)] states that

[TABLE]

holds for all $C^{1}$-vector fields ${\boldsymbol{V}}$. We split ${\boldsymbol{V}}$ into its normal and tangential components according to ${\boldsymbol{V}}=({\boldsymbol{V}}^{\top}{\boldsymbol{n}})\,{\boldsymbol{n}}+\sum_{i=1}^{2}({\boldsymbol{V}}^{\top}{\boldsymbol{\xi}}_{i})\,{\boldsymbol{\xi}}_{i}$ and arrive at

[TABLE]

In the last step we used that ${\boldsymbol{V}}$ is normal and thus ${\boldsymbol{V}}^{\top}{\boldsymbol{\xi}}_{i}=0$ holds. ∎

3.2.  Properties of the Total Variation of the Normal

As part of this section, we seek to establish shape differentiability of our novel objective (1.2). To this end, we utilize that (1.2) is a composition of smooth functions except in the presence of flat regions of positive measure on $\Gamma$. Hence, we first present some results of the material derivatives for the quantities involved, which by themselves are also interesting for a wide variety of other problems.

With respect to the outer normal ${\boldsymbol{n}}$, we first note that under an appropriate regularity assumption, the material derivative exists and is given by

[TABLE]

This result can be found, for instance, in [SokolowskiZolesio1992, Eq. (3.168)] or [Schmidt2010, Lemma 3.3.6]. Notice that $\textup{d}{\boldsymbol{n}}[{\boldsymbol{V}}]$ is tangential because

[TABLE]

Theorem 3.2**.**

Suppose that the orthonormal basis components ${\boldsymbol{\xi}}_{1}$ and ${\boldsymbol{\xi}}_{2}$ are smooth vector fields on a relatively open part $\Gamma_{0}\subset\Gamma$. Then they are shape differentiable and their material derivatives are given by

[TABLE]

The requirement that ${\boldsymbol{\xi}}_{1,2}$ be smooth is not restrictive. The asymmetry in (3.11) stems from the fact that we chose to transform the first basis vector ${\boldsymbol{\xi}}_{1}$ along with (3.1) and then to orthogonalize the second basis vector ${\boldsymbol{\xi}}_{2}$ w.r.t. the first.

Proof.

The proof is by construction. Beginning from a local parametrization ${\boldsymbol{h}}$ of the surface, we give an explicit formula for the tangent basis. The perturbed surface $\Gamma_{\varepsilon}$ is then expressed via a perturbed parametrization ${\boldsymbol{h}}_{\varepsilon}\coloneqq{\boldsymbol{T}}_{\varepsilon}\circ{\boldsymbol{h}}$, where ${\boldsymbol{T}}_{\varepsilon}$ is given by (3.1). 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. $\varepsilon=0$.

Let ${\boldsymbol{h}}:U\subset{}\mathbb{R}^{2}\to\mathbb{R}^{3}$ be a local smooth orthogonal parametrization of $\Gamma_{0}$, i.e., the derivative $D{\boldsymbol{h}}$ is a matrix with orthonormal columns, such that ${\boldsymbol{s}}\in\Gamma$ is locally given by ${\boldsymbol{s}}={\boldsymbol{h}}({\boldsymbol{x}})$ for some ${\boldsymbol{x}}\in U$. Hence, we can define a smooth, orthonormal set of tangent vectors ${\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2}$ via

[TABLE]

where ${\boldsymbol{e}}_{i}$ is the $i$-th canonical basis vector of $\mathbb{R}^{3}$. With respect to ${\boldsymbol{\xi}}_{1}$, we arrive at the normalized tangent vector of the perturbed surface as

[TABLE]

Regarding ${\boldsymbol{\xi}}_{2}$, we proceed in a similar way, but have to apply a Gram–Schmidt step to obtain an orthonormal set of perturbed tangent vectors. Hence, ${\boldsymbol{\xi}}_{2,\varepsilon}$ is given by

[TABLE]

A straightforward differentiation with respect to $\varepsilon=0$ results in the material derivatives given in (3.11). ∎

Theorem 3.3**.**

Under the assumptions of the previous theorem, the derivative $D_{\Gamma}{\boldsymbol{n}}$ of the normal is shape differentiable. The material derivatives of the directional derivatives of ${\boldsymbol{n}}$ in the directions of ${\boldsymbol{\xi}}_{1,2}$ are

[TABLE]

Proof.

With the material derivative of both the normal (3.9) and tangent (3.11) at hand, we apply the chain rule and due to the relationship between spatial and material derivatives (3.6) we arrive at

[TABLE]

Because $\textup{d}{\boldsymbol{\xi}}_{i}[{\boldsymbol{V}}]$ is tangential due to (3.11), the reduction of the full derivative $D$ to the intrinsic derivative $D_{\Gamma}$ on the surface is straightforward provided that ${\boldsymbol{n}}$ is assumed to be extended constantly into a tubular neighborhood of $\Gamma$, i.e., $(D{\boldsymbol{n}}){\boldsymbol{n}}={\boldsymbol{0}}$. ∎

Theorem 3.4**.**

Suppose that the principal curvatures do not vanish simultaneously on $\Gamma$ except possibly on a set of measure zero. Then (1.2) is shape differentiable.

Proof.

Using [SokolowskiZolesio1992, Sec 2.18 and Sec 2.33], shape differentiability of (1.2) can be established for $\Gamma$ of class $C^{1}$ if the integrand of (1.2), i.e.,

[TABLE]

fulfills $g(\varepsilon,\cdot)\in L^{1}(\Gamma_{\varepsilon})$. Likewise, the material and local derivatives have to satisfy $\textup{d}g[{\boldsymbol{V}}]\in L^{1}(\Gamma)$ and $g^{\prime}[{\boldsymbol{V}}]\in L^{1}(\Gamma)$ for all sufficiently smooth vector fields ${\boldsymbol{V}}$ with compact support in the hold-all $\mathcal{D}$.

Since we consider $\Gamma$ to be a smooth surface and ${\boldsymbol{V}}\colon\mathcal{D}\rightarrow\mathbb{R}^{3}$ a smooth vector field, $\Gamma_{\varepsilon}\coloneqq T_{\varepsilon}[{\boldsymbol{V}}](\Gamma)$ is smooth. Moreover, since $\Gamma_{\epsilon}$ is compact, $g(\varepsilon,\cdot)$ is bounded, and one easily deduces $g(\varepsilon,\cdot)\in L^{1}(\Gamma_{\varepsilon})$.

The shape differentiability of the tangent basis is considered in Theorem 3.2 and of the derivative of the normal in Theorem 3.3. Hence, we establish the material derivative of the expression under the square root in (3.16) via the chain rule for a composition of smooth functions. Notice that the local character of the results in Theorems 3.2 and 3.3 is sufficient since (3.16) is independent of the choice of the orthonormal basis, as established in Remark 2.1. Since, by assumption, both principal curvatures do not vanish simultaneously, we have $g\neq 0$ almost everywhere and we arrive at

[TABLE]

from where we conclude $\textup{d}g[{\boldsymbol{V}}]\in L^{1}(\Gamma)$. To establish $g^{\prime}[{\boldsymbol{V}}]\in L^{1}(\Gamma)$, we utilize (3.4). To this end, we extend $g$ constantly in normal direction into a tubular neighborhood of $\Gamma$. As a composition of smooth functions, we conclude $g^{\prime}[{\boldsymbol{V}}]\in L^{1}(\Gamma)$. ∎

Remark 3.5**.**

As per [SokolowskiZolesio1992, Eq. (2.172)], the requirement $g^{\prime}[{\boldsymbol{V}}]\in L^{1}(\Gamma)$ can be omitted if one is only interested in the representation (3.3) of the shape derivative and not in formulation (3.5).

We are now in the position to address the minimization of (1.2). In view of lemma 2.3, this is meaningful only when additional terms are present which prevent the degeneration of the surface to a point. We choose to impose a constraint on the surface area here. We have the following partial result.

Theorem 3.6**.**

Spheres are stationary points for (1.2) among all surfaces $\Gamma$ of constant area.

Proof.

We consider the minimization of (1.2) or equivalently, (2.2), subject to the constraint that the surface area equals the constant $A_{0}>0$. The Lagrangian associated with this problem is given by

[TABLE]

Here $\mu\in\mathbb{R}$ is the Lagrange multiplier to be determined below. The differentiability of the first summand has been considered in Theorem 3.4. On the perturbed domain with surface $\Gamma_{\varepsilon}$ with the perturbation according to (3.1), the Lagrangian reads

[TABLE]

We use the same abbreviation as before in (3.16). The above integral is of type (3.2) and its shape derivative at the unperturbed surface, according to (3.5), is given by

[TABLE]

because $\mu$ is a constant.

When $\Omega$ is a sphere of radius $r$, we are going to show that $\textup{d}\mathcal{L}(0,\mu)[{\boldsymbol{V}}]=0$ holds for all perturbation fields ${\boldsymbol{V}}$ in normal direction and with $\mu=-1/(\sqrt{2}\,r)$. In this setting, the principal curvatures are $k_{1}({\boldsymbol{s}})=k_{2}({\boldsymbol{s}})\equiv 1/r$; see for instance [GrayAbbenaSalamon2006, Chapter 13]. Consequently, g({\boldsymbol{s}})=\mathopen{}\mathclose{{}\left(k_{1}^{2}({\boldsymbol{s}})+k_{2}^{2}({\boldsymbol{s}})}\right)^{\frac{1}{2}}\equiv\sqrt{2}/r is spatially constant and thus $Dg\equiv 0$ holds. We obtain from (3.5)

[TABLE]

Hence, by (3.4), we also have $\textup{d}g[{\boldsymbol{V}}]=g^{\prime}[{\boldsymbol{V}}]$. Using (3.17), we arrive at

[TABLE]

In order to complete the proof of Theorem 3.6, we need to show that

[TABLE]

holds with $c_{0}=-\frac{\sqrt{2}}{r^{2}}$. To this end, we need a tangential Stokes formula as given in Lemma 3.1.

We shall also utilize that $g({\boldsymbol{s}})=\sqrt{2}/r$ is a constant on the sphere of radius $r$. Finally, we utilize

[TABLE]

and thus $(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}={\boldsymbol{\xi}}/r$ holds for $i=1,2$.

The three terms contributing to the material derivative $\textup{d}[(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{i}][{\boldsymbol{V}}]$ in (3.19) are given in (3.15) and we consider them individually. We utilize that the Riemannian metric on $\SS^{2}$ is the Euclidean inner product of the ambient $\mathbb{R}^{3}$, i.e., $\mathfrak{g}({\boldsymbol{a}},\,{\boldsymbol{b}})={\boldsymbol{a}}^{\top}{\boldsymbol{b}}$.

First Term. The insertion of the first term in (3.15) into the left hand side of (3.19) leads to the expression

[TABLE]

The last step follows from (3.8) with $c=1$. Recall from (3.10) that $\textup{d}{\boldsymbol{n}}[{\boldsymbol{V}}]$ is tangential.

Second Term. Inserting the second term in (3.15) into the left hand side of (3.19) leads to the expression

[TABLE]

Third Term. Finally, inserting the third term in (3.15) into the left hand side of (3.19) yields

[TABLE]

The first summand ($i=1$) leads to

[TABLE]

For the second summand ($i=2$), we get one additional term:

[TABLE]

Hence expression (3.23) is zero. Collecting terms (3.21)–(3.23), we have shown that the left hand side in (3.19) amounts to

[TABLE]

Therefore, (3.19) is fulfilled with

[TABLE]

As a consequence of (3.19) we obtain the representation

[TABLE]

for all perturbation fields ${\boldsymbol{V}}$ parallel to ${\boldsymbol{n}}$. Clearly, the term in brackets vanishes when $\mu=-1/(\sqrt{2}\,r)$ holds. This concludes the proof. ∎

Remark 3.7**.**

1. (1)

A numerical study shows that among all ellipsoids of equal area, the sphere is indeed the unique minimizer of (1.3). 2. (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 $\mu=-\sqrt{2}/r^{2}$.

3.3.  The Case of Curves

When $\Omega\subset\mathbb{R}^{2}$ and $\Gamma$ is a one-dimensional manifold, (1.2) and thus (2.2) reduce to the total absolute curvature $\int_{\Gamma}\lvert k\rvert\,\textup{d}{\boldsymbol{s}}$, where $k$ is the single curvature. It is well known that this integral has a minimal value of $2\pi$, which is attained precisely for the boundaries $\Gamma$ of convex, smoothly bounded domains $\Omega\subset\mathbb{R}^{2}$; see [Chen2000, Chapter 21.1] or [BrookBrucksteinKimmer2005]. This case is thus different in two aspects from our setting $\Omega\subset\mathbb{R}^{3}$. On the one hand, $\int_{\Gamma}\lvert k\rvert\,\textup{d}{\boldsymbol{s}}$ is invariant to scale while (1.2) is not, as was shown in (2.5). On the other hand, (1.2) appears to have a much smaller set of minimizers; see Theorem 3.6.

4   Split Bregman Iteration

In this section we propose an Alternating Direction Method of Multipliers (ADMM) iteration, which generalizes the split Bregman algorithm for total variation problems proposed in [GoldsteinOsher2009]. As is well known, ADMM methods introduce a splitting of variables so that minimization over individual variables becomes efficient.

There is very little prior work on ADMM involving manifolds. We are aware of [LaiOsher2014, KovnatskyGlashoffBronstein2016, ZhangMaZhang2017_preprint, WangYinZeng2018], all of which utilize particular manifolds and their embeddings into some vector space in order to formulate the splitting constraint. By contrast, in our setting the constraint is formulated pointwise in the tangent bundle of $\SS^{2}$. In addition, and even though we do not emphasize this aspect throughout the paper, the primary variable $\Omega$ in problem (4.1) lives on a manifold of shapes. That said, we will not attempt a convergence proof for the proposed split Bregman iteration here but leave it for future research.

In our setting, the primary variable is the unknown domain $\Omega$. Notice that $\Omega$ also determines its boundary $\Gamma$ as well as its normal vector field ${\boldsymbol{n}}$, and we will always consider $\Gamma$ and ${\boldsymbol{n}}$ as a function of $\Omega$. The splitting is achieved through the introduction of a new variable ${\boldsymbol{d}}$, which is independent of $\Omega$, $\Gamma$ and ${\boldsymbol{n}}$. At the solution, we require the coupling condition ${\boldsymbol{d}}=D_{\Gamma}{\boldsymbol{n}}$ to hold across $\Gamma$. We recall that $D_{\Gamma}{\boldsymbol{n}}$ denotes the derivative (push-forward) of ${\boldsymbol{n}}$. At the point ${\boldsymbol{s}}\in\Gamma$, $(D_{\Gamma}{\boldsymbol{n}})({\boldsymbol{s}})$ maps $\mathcal{T}_{{\boldsymbol{s}}}\Gamma$ into $\mathcal{T}_{{\boldsymbol{n}}({\boldsymbol{s}})}\SS^{2}$.

Written in terms of $\Omega$ and the secondary variable ${\boldsymbol{d}}=({\boldsymbol{d}}_{1},{\boldsymbol{d}}_{2})\colon\Gamma\to\mathcal{T}_{{\boldsymbol{n}}(\cdot)}\SS^{2}$, problem (1.3) becomes

[TABLE]

Notice that for convenience of notation, we represent $D_{\Gamma}{\boldsymbol{n}}$ in terms of its actions on the two basis vectors ${\boldsymbol{\xi}}_{i}$.

Note also that at the point ${\boldsymbol{s}}\in\Gamma$, the equality ${\boldsymbol{d}}_{i}=(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{i}$ is in the tangent space $\mathcal{T}_{{\boldsymbol{n}}({\boldsymbol{s}})}\SS^{2}$. We therefore introduce Lagrange multipliers ${\boldsymbol{\lambda}}=({\boldsymbol{\lambda}}_{1},{\boldsymbol{\lambda}}_{2})$, belonging to the same space, and define the augmented Lagrangian associated with (1.3) as follows,

[TABLE]

After the usual re-scaling ${\boldsymbol{b}}_{i}\coloneqq{\boldsymbol{\lambda}}_{i}/\lambda$, we can rewrite (4.2) as

[TABLE]

The main difference to an ADMM method in Euclidean or Hilbert spaces is that the vector fields ${\boldsymbol{d}}_{i}$ and ${\boldsymbol{b}}_{i}$ have values in the tangent space $\mathcal{T}_{{\boldsymbol{n}}(\cdot)}\SS^{2}$. Hence they must be updated whenever $\Omega$ and thus the normal vector field ${\boldsymbol{n}}$ are changing.

We outline our proposed method for (4.1) as Algorithm 4.1. As expected for methods of the ADMM class, we successively optimize with respect to the variables $\Omega$ and $({\boldsymbol{d}}_{1},{\boldsymbol{d}}_{2})$ independently and then perform a simple update step for the multiplier $({\boldsymbol{b}}_{1},{\boldsymbol{b}}_{2})$. We address each of these steps in the following subsections.

Algorithm 4.1**.**

Split Bregman method for (4.1)

0:*  Initial domain $\Omega^{(0)}$ *

0:*  Approximate solution of (4.1) *

1:*  Set ${\boldsymbol{b}}^{(0)}\coloneqq{\boldsymbol{0}}$, ${\boldsymbol{d}}^{(0)}\coloneqq{\boldsymbol{0}}$ *

2:*  Set $k\coloneqq 0$ *

3:*  while not converged do

4:*     Perform some gradient steps for $\Omega\mapsto\mathcal{L}(\Omega,{\boldsymbol{d}}^{(k)},{\boldsymbol{b}}^{(k)})$ starting from $\Omega^{(k)}$, to obtain $\Omega^{(k+1)}$ *

5:*     Parallely transport the multiplier estimate ${\boldsymbol{b}}^{(k)}$ pointwise from $\mathcal{T}_{{\boldsymbol{n}}^{(k)}(\cdot)}\SS^{2}$ to $\mathcal{T}_{{\boldsymbol{n}}^{(k+1)}(\cdot)}\SS^{2}$ along the geodesic from ${\boldsymbol{n}}^{(k)}$ to ${\boldsymbol{n}}^{(k+1)}$ *

6:*     Parallely transport the basis vectors ${\boldsymbol{\xi}}_{i}$ pointwise from $\mathcal{T}_{{\boldsymbol{n}}^{(k)}(\cdot)}\SS^{2}$ to $\mathcal{T}_{{\boldsymbol{n}}^{(k+1)}(\cdot)}\SS^{2}$ along the geodesic from ${\boldsymbol{n}}^{(k)}$ to ${\boldsymbol{n}}^{(k+1)}$ for $i=1,2$ *

7:*     Set ${\boldsymbol{d}}^{(k+1)}\coloneqq\operatorname*{arg\,min}\mathcal{L}(\Omega^{(k+1)},{\boldsymbol{d}}^{(k)},{\boldsymbol{b}}^{(k)})$, see (4.7) *

8:*     Update the Lagrange multipliers, i.e., set ${\boldsymbol{b}}_{i}^{(k+1)}\coloneqq{\boldsymbol{b}}_{i}^{(k)}+(D_{\Gamma}{\boldsymbol{n}}^{(k+1)})\,{\boldsymbol{\xi}}_{i}-{\boldsymbol{d}}_{i}^{(k+1)}$ for $i=1,2$ *

9:*     Set $k\coloneqq k+1$ *

10:*  end *while

4.1.  The Shape Optimization Step

We first address the minimization of (4.3) w.r.t. the shape ${\Omega}$, while the quantities ${\boldsymbol{d}}_{1},{\boldsymbol{d}}_{2},{\boldsymbol{b}}_{1},{\boldsymbol{b}}_{2}$ are fixed, or, more precisely, passively transformed along with $\Omega$. The main effort is to calculate the shape derivative of (4.3).

The derivative of the first term in (4.3), i.e., $\textup{d}\ell(u({\Omega}),{\Omega})[{\boldsymbol{V}}]$, is not specified here since it depends on the particular PDE underlying the solution operator $u(\Omega)$ and the loss function $\ell$ considered. This derivative can be obtained by standard shape calculus techniques, which are not our concern here. A concrete example is considered in the companion paper [BergmannHerrmannHerzogSchmidtVidalNunez2019:2_preprint].

Next we consider the second term in (4.3). Due to the chosen splitting, the vector fields ${\boldsymbol{d}}_{i}$ are merely transformed along with $\Omega$ according to the perturbation (3.1) and thus we define their perturbed counterparts as

[TABLE]

and likewise for ${\boldsymbol{b}}_{i}$ and \mathopen{\bigl{(}}{\lvert{\boldsymbol{d}}_{1}\rvert_{\mathfrak{g}}^{2}+\lvert{\boldsymbol{d}}_{2}\rvert_{\mathfrak{g}}^{2}}\mathclose{\bigr{)}}^{1/2}. As a consequence, their material derivatives vanish and the directional derivative of the second term of (4.3) becomes

[TABLE]

Finally we address the terms \displaystyle\int_{{\Gamma}}\bigl{\lvert}{\boldsymbol{d}}_{i}-({D_{\Gamma}{\boldsymbol{n}}})\,{{\boldsymbol{\xi}}_{i}}-{\boldsymbol{b}}_{i}\bigr{\rvert}_{\mathfrak{g}}^{2}\,{\textup{d}{\boldsymbol{s}}}, $i=1,2$. The vector fields ${\boldsymbol{d}}_{i}$ and ${\boldsymbol{b}}_{i}$ are transformed according to (4.4) and thus we need not consider their material derivatives. However, we do need to track the dependencies of $(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{i}$. The respective shape derivative is given in (3.15).

We summarize these findings in the following theorem.

Theorem 4.2**.**

Suppose that the $({\boldsymbol{d}}_{1},{\boldsymbol{d}}_{2})$ do not vanish simultaneously on $\Gamma$ except possibly on a set of measure zero, and that the loss term $\ell(u(\Omega),\Omega)$ is shape differentiable. Then the augmented Lagrangian (4.3) is shape differentiable and its shape derivative is given by

[TABLE]

with \textup{d}\mathopen{\bigl{(}}{(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{i}}\mathclose{\bigr{)}}[{\boldsymbol{V}}] from (3.15).

The shape derivative in (4.2) is the basis of any shape optimization procedure. After all, the minimization of (4.3) w.r.t. the domain $\Omega$ represents a fairly standard shape optimization problem. We convert the shape derivative (4.2) into a shape gradient vector field ${\boldsymbol{U}}$ by means of an appropriate inner product. Then, instead of minimizing (4.3) w.r.t. $\Omega$ to a certain accuracy, in practice we only perform a few gradient steps per ADMM iteration. This is in line with [GoldsteinOsher2009], where a Gauss–Seidel sweep is proposed instead of an exact solve.

Still, the terms in (3.15) would be tedious to implement by hand. In our implementation, which is detailed in the companion paper [BergmannHerrmannHerzogSchmidtVidalNunez2019:2_preprint], the shape derivative (4.2) is conveniently evaluated by algorithmic differentiation techniques on the discrete level.

4.2.  The Total Variation Minimization Step

Before addressing the minimization of (4.3) w.r.t. ${\boldsymbol{d}}$ we must note that the data ${\boldsymbol{b}}_{i}$ at any point ${\boldsymbol{s}}\in\Gamma$ has to belong to the tangent space $\mathcal{T}_{{\boldsymbol{n}}({\boldsymbol{s}})}\SS^{2}$. Since the surface $\Gamma$ and hence the field of normal vectors is changing during the shape optimization step, we must first move the data ${\boldsymbol{b}}_{i}$ into the new tangent space. This is achieved via parallel transport along geodesics on $\SS^{2}$. Suppose for brevity of notation that ${\boldsymbol{n}}^{-}$ denotes the old normal vector field along the boundary $\Gamma^{-}$ of the previous iterate $\Omega^{-}$. Then ${\boldsymbol{b}}_{i}^{-}\in\mathcal{T}_{{\boldsymbol{n}}^{-}(\cdot)}\SS^{2}$ needs to be parallely transported into ${\boldsymbol{b}}_{i}\in\mathcal{T}_{{\boldsymbol{n}}(\cdot)}\SS^{2}$ along the geodesic from ${\boldsymbol{n}}^{-}$ to ${\boldsymbol{n}}$. This step is inexpensive since geodesics and parallel transport on $\SS^{2}$ are available in terms of explicit formulas; see Appendix A. Since we explicitly refer to them, also the basis vectors ${\boldsymbol{\xi}}_{i}^{-}$ need to be parallely transported into ${\boldsymbol{\xi}}_{i}$ in the same way as above.

We can now address the minimization of (4.3) w.r.t. the field ${\boldsymbol{d}}=({{\boldsymbol{d}}_{1}},{{\boldsymbol{d}}_{2}})$. Since the first term in (4.3) does not depend on ${\boldsymbol{d}}$, we are left with the minimization of

[TABLE]

where ${{\boldsymbol{d}}_{1}},{{\boldsymbol{d}}_{2}}$ are sought pointwise in the tangent spaces $\mathcal{T}_{{\boldsymbol{n}}(\cdot)}\SS^{2}$. We recall that the latter are two-dimensional subspaces of $\mathbb{R}^{3}$. At this point it is important to note that the data $(D_{\Gamma}{\boldsymbol{n}})\,{\boldsymbol{\xi}}_{i}+{\boldsymbol{b}}_{i}$ belongs to the same tangent spaces. Therefore, just like in the Euclidean setting, the minimizer of (4.6) can be evaluated explicitly and inexpensively via a pointwise, vectorial shrinkage operation, i.e.,

[TABLE]

Here we abbreviated

[TABLE]

4.3.  The Multiplier Update

An update of the Lagrange multiplier fields $({{\boldsymbol{b}}_{1}},{{\boldsymbol{b}}_{2}})$ is achieved, analogously to the Euclidean setting, by replacing ${{\boldsymbol{b}}_{i}}$ with

[TABLE]

Notice again that all quantities belong to the subspace $\mathcal{T}_{{\boldsymbol{n}}(\cdot)}\SS^{2}$ of $\mathbb{R}^{3}$.

5   Conclusions and Outlook on the Discrete Setting

In this paper we introduced an analogue of the total variation prior for the normal vector field (1.2) defined on the boundary $\Gamma$ of smooth domains $\Omega\subset\mathbb{R}^{3}$. This functional is also known as the total root mean square curvature (2.2). We have shown in theorem 3.6 that it admits spheres as stationary points under an area constraint and we conjecture that spheres are in fact global minimizers.

We proposed a split Bregman (ADMM) scheme for the numerical solution of shape optimization problems (1.3) involving the total variation of the normal. In contrast to a Euclidean ADMM, as proposed for instance in [GoldsteinOsher2009], the normal vector data belongs to the sphere $\SS^{2}$. Therefore, the formulation of the ADMM method requires concepts from differential geometry. An analysis of the Riemannian ADMM scheme is beyond the scope of this paper and will be presented elsewhere.

In the companion paper [BergmannHerrmannHerzogSchmidtVidalNunez2019:2_preprint], we consider a discrete version of the total variation of the normal functional (1.2), which applies to piecewise flat surfaces. These arise naturally when surfaces are represented by meshes, or when they are discretized for computational purposes. We also discuss the utility of the discrete TV of the normal as a shape prior to recover piecewise flat surfaces in geometric inverse problems, including PDE constraints.

Appendix A The Sphere as a Riemannian Manifold

In this section we provide some useful formulas for the sphere

[TABLE]

equipped with the Riemannian metric obtained from the pull back of the Euclidean metric from the ambient space $\mathbb{R}^{3}$. We are going to represent points ${\boldsymbol{n}}\in\SS^{2}$ by vectors in $\mathbb{R}^{3}$. Moreover, we identify the tangent space at ${\boldsymbol{n}}$ with the two-dimensional subspace

[TABLE]

We utilize the Riemannian metric $\mathfrak{g}({\boldsymbol{a}},\,{\boldsymbol{b}})={\boldsymbol{a}}^{\top}{\boldsymbol{b}}$ in $\mathcal{T}_{{\boldsymbol{n}}}\SS^{2}$ and the norm $\lvert{\boldsymbol{a}}\rvert_{\mathfrak{g}}=({\boldsymbol{a}}^{\top}{\boldsymbol{a}})^{1/2}$.

The geodesic distance between any two ${\boldsymbol{n}},{\boldsymbol{n}}^{\prime}\in\SS^{2}$ is given by

[TABLE]

The geodesic curve $\gamma(\,\cdot\,;{\boldsymbol{n}},{\boldsymbol{\xi}})\colon\mathbb{R}\to\SS^{2}$ departing from ${\boldsymbol{n}}\in\SS^{2}$ in the direction of ${\boldsymbol{\xi}}\in\mathcal{T}_{{\boldsymbol{n}}}\SS^{2}$ is given by

[TABLE]

The exponential map is thus given by

[TABLE]

The logarithmic map is the inverse of the exponential map w.r.t. to the tangent direction ${\boldsymbol{\xi}}$. In other words, ${\boldsymbol{\xi}}=\log_{{\boldsymbol{n}}}{\boldsymbol{n}}^{\prime}$ holds if and only if ${\boldsymbol{\xi}}$ is the unique element in $\mathcal{T}_{{\boldsymbol{n}}}\SS^{2}$ such that $\exp_{{\boldsymbol{n}}}{\boldsymbol{\xi}}={\boldsymbol{n}}^{\prime}$. The logarithmic map is well-defined whenever ${\boldsymbol{n}}\neq-{\boldsymbol{n}}^{\prime}$ holds. In this case, we have

[TABLE]

Finally we require the concept of parallel transport of a tangent vector from one tangent space to another, along the unique shortest geodesic connecting the base points. Specifically, the parallel transport $P_{{\boldsymbol{n}}\to{\boldsymbol{n}}^{\prime}}\colon\mathcal{T}_{{\boldsymbol{n}}}\SS^{2}\to\mathcal{T}_{{\boldsymbol{n}}^{\prime}}\SS^{2}$ along the unique shortest geodesic $\gamma(\,\cdot\,;{\boldsymbol{n}},\log_{{\boldsymbol{n}}}{\boldsymbol{n}}^{\prime})$ connecting ${\boldsymbol{n}}$ and ${\boldsymbol{n}}^{\prime}\neq-{\boldsymbol{n}}$ is given by

[TABLE]

see for instance [HosseiniUschmajew2017] and [Persch2018, Section 2.3.1], repectively. Here we used the abbreviations ${\boldsymbol{v}}=\log_{{\boldsymbol{n}}}{\boldsymbol{n}}^{\prime}$, $\lvert{\boldsymbol{v}}\rvert_{\mathfrak{g}}=d({\boldsymbol{n}},{\boldsymbol{n}}^{\prime})$ and ${\boldsymbol{u}}=\frac{{\boldsymbol{v}}}{\lvert{\boldsymbol{v}}\rvert_{\mathfrak{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

[TABLE]

which holds true since the norm of the logarithmic map is

[TABLE]

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

[TABLE]

Acknowledgments

The authors would like to thank two anonymous reviewers for their constructive criticism which helped improve the paper.

This work was supported by DFG grants HE 6077/10–1 and SCHM 3248/2–1 within the Priority Program SPP 1962 (Non-smooth and Complementarity-based Distributed Parameter Systems: Simulation and Hierarchical Optimization), which is gratefully acknowledged.