Total Roto-Translational Variation

We consider curvature depending variational models for image regularization, such as Euler's elastica. These models are known to provide strong priors for the continuity of edges and hence have important applications in shape-and image processing. We consider a lifted convex representation of these models in the roto-translation space: In this space, curvature depending variational energies are represented by means of a convex functional defined on divergence free vector fields. The line energies are then easily extended to any scalar function. It yields a natural generalization of the total variation to the roto-translation space. As our main result, we show that the proposed convex representation is tight for characteristic functions of smooth shapes. We also discuss cases where this representation fails. For numerical solution, we propose a staggered grid discretization based on an averaged Raviart-Thomas finite elements approximation. This discretization is consistent, up to minor details, with the underlying continuous model. The resulting non-smooth convex optimization problem is solved using a first-order primal-dual algorithm. We illustrate the results of our numerical algorithm on various problems from shape-and image processing.


Introduction
It was observed at least since [58] that line energies such as Euler's "Elastica" (here for a smooth curve γ): γ [1+] κ 2 ds could be natural regularizers for the completion of missing contours in images. This idea was based in particular on observations of G. Kanizsa [46,47] about the way our perception can "invent" apparent contours. The model in [58] (see also [40,39] for interesting attempts to solve it with phase-field methods) was variational in nature and since then, many attempts have been proposed to study and address the minimization of such line energies, both theoretically and numerically.
From the theoretical point of view, the lower-semicontinuity of these energies already is a challenge, which has been studied in many papers (in connexion to the applications to computer Riemannian and sub-Finslerian geodesic tracking in the roto-translation group [11,34]. As said, we wish to extend these variational methods to level sets representations. The family of problems we are interested in is introduced in a paper by Masnou and Morel [52] (see also [53]), which addresses the problem of image inpainting. Their initial idea is to minimize, in the inpainting domain D, an energy of the form D

+ div
Du |Du| p |Du| for some p ≥ 1, where u is a bounded variation function and appropriate boundary conditions are given. A theoretical study of this energy, which in general is not lower-semicontinuous (lsc), is found in [4], in particular it is shown that if p > 1 (in dimension 2), it is lsc on C 2 functions. Interestingly, also this study relies on a space/direction representation (and more precisely on varifolds). An interesting co-area formula for the relaxed envelope is also shown in [51].
There have been many attempts to numerically solve Masnou and Morel's model, which is very difficult to tackle, being highly non convex. Most of the authors introduce auxiliary variables [9,26,10], for instance representing the orientation, which is already close to the idea of the Gauss map or Roto-translational representation. Recent techniques based on Augmented Lagrangian methods (for coupling the auxiliary variables) [74,80,78,79,43,42,8] have shown to be quite efficient, despite the lack of convexity and hence convergence guaranties.
What we propose here is to rather introduce a functional which may be sees as a convex relaxation of Masnou and Morel's. The method we propose is based on two ingredients: A lifting in the Roto-translation space of curves which allows to write the Elastica energy (or any convex function of the curvature) as a convex function, as in classical Gauss map based approaches, and (formally) a decomposition of functions as sum of characteristic functions of sets with finite energy. This allows to define a convex functional which is defined on grey-level functions and penalises the curvature of the level lines. It is however easy to check that it is in general strictly below Masnou and Morel's functional, in particular a function of finite energy needs not have its level sets regular (in practice, they could be intersections of regular sets, in addition to having possible cusps, for the same reasons as in [12]), see for example Figures 5 and 9(e). The only theoretical result which we can show is that C 2 curves are tightly represented by our convex relaxation.
There is a close relationship between our approach and the functionals in [22] ("TVX") and [23], based on similar representations (but [22] omits to preserve the boundary of the lifted current). In fact, the functional we build is, as we show further on, a new expression of the previous convex relaxation of the Elastica energy (and variants) proposed in [23]. This stems from the identity of the corresponding dual problems. Although it boils down to the same energy, it is introduced in a much simpler way, as the primal expression in [23] requires to work in a space where the point, the tangent and the curvature are lifted as independent variables. As a consequence, we also can provide a simpler discretization, and the tightness of the relaxation for C 2 sets was unnoticed in [23]. Eventually, we should mention that part of what we propose here could be generalized to arbitrary dimension or co-dimension, following the techniques in [6]. However, computationally, our construction for 1D curves in the plane already needs to work in a 3D space, and more complex models seem at this moment intractable. Our future research will rather focus on improving the discretization and the optimization of the bidimensional case.
The paper is organised as follows: In Section 2 we describe the lifting and introduce our functional. We state our main result, which shows that characteristic functions of C 2 sets are well represented by our convex relaxation (Theorem 1). In the next Section 3 we show an approximation result with smooth functions, and compute a dual representation of the functional (which shows it is identical to [23]). Then, in Section 4, we describe how the functional is discretized and show a few experimental results. The last Section 5 is devoted to the proof if Theorem 1. Some technical tools are found in the Appendix, in particular, Appendix A shows, up to (hopefully minor) transformations, the consistency of our implementation in Section 4.

The functional
We consider Ω a bounded open domain in the plane, and E ⊂ Ω a set with C 2 (and to simplify, connected) boundary. Assume that this boundary ∂E is given by a parameterized curve (x 1 (t), x 2 (t)) with parameter t ∈ [0, 1]. The (extrinsic) curvature κ E of ∂E is classically defined as the ratio between the variation of the tangential angle θ and the variation of its arc length s, that is The main idea of the lifting is now to consider a higher dimensional representation of the parametric curve in the 3D roto-translation (RT) space, which is obtained by adding the tangential angle θ as an additional dimension to Ω. In this space, we now consider a parameterized 3D curve (x 1 (t), x 2 (t), θ(t)) which lifts the boundary ∂E to the roto-translational space. Figure 1 shows an example where we lift the boundary of a disk (which is a circle) to the RT space. Observe that the lifted boundary is represented by a 3D helix. Now, we define for all t ∈ [0, 1] the tangential vector p(t) = (p x (t), p θ (t)) with The curvature is therefore given by In this work we consider f : R → [0, +∞] a convex, lsc function and want to define a convex lsc extension (to grey-level valued functions) of energies of the type where E ⊂ Ω is a set with C 2 boundary, and κ E is the curvature of the set. Using our tangential vector p(t), it is easy to see that the energy can be (formally) written as 1]) is the lifted curve and τ (x, θ) its normalized tangential vector, given by τ (x(t), θ(t)) = p(t)/|p(t)| for all t ∈ [0, 1]. A precise definition of the energy will be given below.
Let us briefly discuss three instances of energies that will typically appear in applications. In all cases α > 0 will be a tuning parameter that can be used to balance the influence of the curvature term with respect to the length. Figure 1: The gray line is the lifting of the boundary ∂E of the disk E to the roto-translational space Ω × S 1 .
1. f 1 (t) = 1 + α|t|. This energy penalizes the arclength plus the absolute curvature, hence we might expect that this type of energy will allow also for corners. This type of energy has been studied with a different approach in [22]. The energy is given by: 2. f 2 (t) = 1 + α 2 |t| 2 . This energy penalizes the arclength of the lifted curve in the RT space which, in some sense, corresponds to the "Total Roto-translational Variation". It yields the energy: which is nothing but the length of the lifted curve Γ E in a Riemanian metric. This functional was considered for perceptional completion problems probably first in [65]. It coincides with sub-Riemannian [65,33,20,11] and more precisely the sub-Finslerian models in [34] (with the constraint of positive direction on the velocity).
Observe that the second term in the energy is a "quadratic over linear" function which is still convex, provided τ x is constrained in a half space or a half line (which of course will depend on the θ variable).
In all these examples, we find that the energy which we are interested in ends up represented possibly as a convex function of the measure σ = τ H 1 Γ E . Moreover, this σ is not an arbitrary measure. In particular, it satisfies two important constraints: first, by construction, it is a circulation and has zero divergence in Ω × S 1 (it can have source terms on ∂Ω × S 1 ): indeed, given any smooth function ψ, ∇ψ · σ = Γ E ∂ τ ψdH 1 vanishes if Γ E is a closed curve, or if it ends on the boundary and ψ has compact support. Second, its marginals in Ω × S 1 , which we can formally denote S 1 σ (and which are also divergence free), coincide with a 90 • -rotation of the measure Dχ E . We will show now how to generalize this construction to arbitrary sets or functions (with bounded variation).
In the whole paper, we assume that there exists γ > 0 such that We also introduce the recession function which is a convex, one-homogeneous function possibly infinite on (−∞, 0) or/and (0, +∞). It is easy to check that it is the support function of dom f * , the domain of the convex conjugate f * of f : Indeed, for y ∈ dom f * , s > 0, as s → ∞; on the other hand, as f * ≥ −γ thanks to (1), for s > 0 one has We let then, for p = (p x , p θ ) ∈ R 3 with p x = 0, Here, θ denotes the unit planar vector (cos θ, sin θ) T (by a slight abuse of notation, we will also denote in this way the vector (cos θ, sin θ, 0) T ∈ R 3 ). It is then classical [64, § 13] and easily follows from (2) that one has h(θ, p) = sup ξ · p : that is, the one-homogeneous function h is the support function of the convex set in the right-hand side of (4). Indeed, first of all, the sup in (4) is ∞ if p x is not λθ, λ ≥ 0. While if it is of this form, then (taking the suprema over the ξ's which satisfy ξ Observe that (1) yields that for all p ∈ R 3 , We now introduce the functional Here x ⊥ = (x 2 , −x 1 ) is a 90 • rotation in the plane. The last condition is understood as follows: In particular in general the fields σ appearing in (6) are free-divergence bounded Radon measures in Ω × S 1 with values in R 3 (we denote M 1 (Ω × S 1 ; R 3 ) the space of such measures) and the proper way to write the integral is rather where σ/|σ| is the Radon-Besicovitch derivative of σ w.r. its total variation |σ|. Notice that for any u and any admissible σ, thanks to (5), because of (7), taking the supremum with respect to the test function ϕ with |ϕ(x)| ≤ 1 everywhere. This shows that F (u) bounds the BV seminorm. We denote by K(θ) the closed convex set of R 3 whose support function is h(θ, ·), as already seen it is given by Remark that (1) implies that f * (s) ≤ − γ 2 − s 2 for |s| ≤ γ, implying in particular that 0 is in the interior of K(θ). We then let (following for instance [21,63]) Then, we claim that for any measure σ ∈ M 1 (Ω × S 1 ; R 3 ) showing in particular that this integral is a lower-semicontinuous function of σ (for the weak- * convergence). To prove this claim (which is standard), first observe that if ϕ ∈ K, then for all ψ − ∈ C 0 c (Ω × S 1 ; R + ) and ψ ∈ C 0 c (Ω × S 1 ), also ϕ − ψ − θ + ψθ ⊥ ∈ K. One deduces easily that the supremum in (8) is infinite if one of the measures (θ · σ x ) − or θ ⊥ · σ x does not vanish. If both vanish, it means that σ = (λθ, σ θ ) for some nonnegative measure λ. Hence the supremum becomes and (8) is then deduced in a standard way (with a Besicovitch covering of Ω × S 1 with respect to the measure |σ| by balls where λ, σ θ are "almost constant", and choosing then for ϕ the "right function" in each ball of the covering-the fact that zero is in the interior of K, which allows to multiply functions in K by a cut-off, is important here). We can deduce that also F is lower semicontinuous: consider a sequence (u n ) with F (u n ) ≤ c < ∞. Then u n is bounded in BV and converges (up to a constant and a subsequence) to some u in L 1 (Ω).
If σ n reaches the value of F (u n ) up to 1/n, one has that Ω×S 1 |σ n | ≤ c < ∞ so that σ n * σ (up to a subsequence), as measures, and obviously div σ = 0. Clearly, also (7) passes to the limit. Hence, by lower semicontinuity, This shows that F defines a convex, lower semicontinuous functional on BV (Ω). Our main theoretical result is the following theorem, which shows that if the argument u is the characteristic function of a smooth enough set, then F (u) coincides, as expected, with a curvature-dependent energy of the boundary of the set {u = 1}. The proof of this result is postponed to Section 5.
Remark 2.1. One could hope that F coincides with the lower semicontinuous envelope of its restriction to C 2 sets, with respect to the L 1 convergence. However, simple examples show that it is not the case. Many examples where it fails are found in [23]. In general, we expect the relaxation F to be strictly below the relaxations found in the literature since [12], based on the (L 1 ) approximation of sets with smooths sets. In particular, in the case of [12,Fig. 4.2], our relaxation will certainly be below twice the energy reached by the pattern in [12, Fig. 1 2. If f has growth one and is even, one expects that the same result holds for piecewise C 2 sets. Remark 2.3. We believe that the proof below could be extended to show that if u ∈ BV (Ω) is a function with C 2 level sets, and such that there exists a continuous function κ(x) which coincides with the curvature of ∂{u > s} for all level s ∈ R, then one should have Hence in that case our functional would coincide with Masnou and Morel's [53]. A typical example of a function u for which such a function κ(x) exists is (in a convex domain Ω, see [54]) a solution of the so-called "Rudin-Osher-Fatemi" functional for some λ > 0, in case f is continuous. Then, also u is continuous and κ(x) is given by λ(u−f ). In that case, in addition, κ(x) is also the curvature of the level sets of u = h(u) for any nondecreasing function h such that h(u) is still in BV (Ω).
3 Some properties of the functional 3.1 Approximation by smooth functions Assume Ω is a bounded convex set. Then for any u ∈ L 1 (Ω) with F (u) < ∞, there exists (u n ) n a sequence of functions with u n ∈ C ∞ (Ω), which converge to u in L 1 (Ω), and such that lim Remark 3.2. We believe that, upon replacing C ∞ (Ω) with C ∞ (Ω), the result should be true in any domain.

Dual representation
We can show the following representation for F , which in particular establishes that F coincides exactly with the (more complicated) relaxation R * * el defined in [23].
Assume ∂Ω is connected. Then the functional (6) can also be represented in the following form If ∂Ω is not connected, then same formula holds, however the functions ϕ should be in C 1 (Ω × S 1 ), with ∇ϕ ∈ C 0 c (Ω × S 1 ).
Corollary 3.4. The functional F (u) coincides with the relaxation of the Elastica energy R * * el (u) defined in [23].
Indeed, formula (12) is the same as (35) in [23,Theorem 13]. The primal formulation (6) we introduce here is quite simpler as the original formulation in [23], though, as the "κ" variable is implicit in our formulation (it comes naturally as the derivative of the orientation) and does not need to be lifted.
Proof of Prop. 3.3. We use the standard perturbation approach to duality [35], in the duality (C 0 0 (Ω × S 1 ; R 3 ), M 1 (Ω × S 1 ; R 3 )), where C 0 0 denotes the L ∞ -closure of the functions with compact support (hence, the functions which vanish on the boundary) and M 1 the totally bounded vector-valued Radon measures. Defining, for p ∈ M 1 (Ω × S 1 ; R 3 ), the function so that F (u) = G(0), we observe that G is lower semicontinuous. Indeed if p n * p and σ n is a minimizer in the definition (13) of G(p n ), as before thanks to (5), σ n + p n is bounded and up to a subsequence, converges to a measure σ + p where σ is admissible in (13). It follows that In particular, we deduce that G = G * * so that It remains to compute G * (η), to show that (12) follows: and we find that Recalling that the sup is over the σ admissible, we find that Now, given any σ with div σ = 0, we can define u such that hence the first supremum in (14) is the same as the sup over all σ with vanishing divergence. Hence, it is zero or +∞, depending on whether η − ψ is a gradient or not. We find eventually G * (η) is finite only when there exists ϕ ∈ C 1 (Ω × S 1 ) with ∇ϕ = 0 on ∂Ω × S 1 , such that η = ψ + ∇ϕ. The proposition follows.

Numerical Experiments
For numerical solution of the proposed model, we need to discretize both the 2D image domain as well as the 3D domain of the roto-translation space. Due to the high anisotropy of the energy, one has to be extremely careful in the choice of the discretization scheme in order to preserve a maximal degree of rotational invariance while keeping the numerical diffusion as low as possible. It turns out that a 2D-3D staggered grid version of an averaged first-order Raviart-Thomas divergence conforming discretization [62] yields the best results. More elaborate discretizations using for example adaptive grids or higher-order approximations will be subject for future study. Figure 2 shows a visualization of our combined 2D-3D staggered grid approach. The first grid is given by a 2D grid of pixels discretizing the image domain Ω. The second grid is given by a 3D The bottom grid represents the 2D image, and the image intensities u are stored at the vertices of the squares S i,j . The gray arrows indicate differences between the pixels. The volume V i,j,k represents one element of the corresponding discretization of the RT space. The black arrows through the faces of the volume represent the 3D vector field σ.

Staggered averaged Raviart-Thomas discretization
grid of volumes which discretize the roto-translation space Ω × S 1 . In the visualization we show a few image pixels of the 2D grid (gray nodes) and one cube of the corresponding 3D grid. First, we start by describing the discretization of the image domain. We assume here that Ω is a square or rectangle, or in general a convex domain, which we will discretize at a scale δ x > 0. We will restrict ourselves to the case of a rectangle to simplify the notation. Then, we let Ω δx denote the discretized image domain which is given by denotes the square centered around the index i = (i, j). The intensity values of the discrete image are stored at the vertices half-grid points (such as the four vertices of S i,j ), and for notational simplicity we introduce the following index set: Hence, the discrete image u ∈ R I consists of N 1 × N 2 discrete pixels with values (u i ) i∈I . We will later also need the following index sets: which index the middle points of the edges between two adjacent pixels along the first and second spatial dimensions. Next, we let δ θ > 0 be the discretization step for the angular variable. We must assume that We then denote by Γ δ with δ = (δ x , δ θ ) the discretized 3D roto-translation space. It is given by denotes the volume centered around the 3D index j = (i, j, k). We associate an index set J of center locations of the volumes V i,j,k : We shall also introduce the following three index sets: which we will use to index the facets of the volumes. Observe that the 3D points which are orthogonal to the first spatial dimension. Similarly, the 3D points (iδ x , (j − 1 2 )δ x , kδ θ ) for all (i, j − 1 2 , k) ∈ J 2 correspond to the middle points of the facets which are orthogonal to the second spatial dimension. Finally, the 3D points (iδ x , jδ x , (k − 1 2 )δ θ ) for all (i, j, k − 1 2 ) ∈ J θ correspond to the middle points of the facets which are orthogonal to the angular dimension. Now, we consider a discrete vector field σ = (σ 1 , respectively. The discrete vector field σ can also be seen as a Raviart-Thomas vector field [62], defined everywhere in Γ δx,δ θ . It is obtained by an affine extension: where ∆(t) = max(0, 1 − |t|) is the usual linear interpolation kernel.
It is well know that such a Raviart-Thomas field has a divergence (defined everywhere in the distributional sense) which is obtained in each volume by summing the finite differences of the average fluxes through each pair of opposite facets. In our model, the discrete field σ is constrained to be divergence free, hence for each volume V j , we must impose that: The roto-translation space is periodic in the θ direction and hence, whenever k = N θ , we identify the index N θ + 1 2 with the index 1 2 in the computation of the last term of the discrete divergence. To simplify our presentation, we introduce a discrete divergence operator D : R J 1 ∪J 1 ∪J θ → R J , so that the divergence free constraint (22) can be compactly written as Next, to implement a discrete version of the consistency condition we observe that averaging (in the direction of θ) the first two components of the field σ yields a rotated (by π/2) gradient which should locally coincide with the finite differences of the discrete image u. This yields the following discrete consistency condition: We introduce a projection operator P : R J 1 ∪J 2 → R I 1 ∪I 2 and a discrete (rotated) gradient operator G : R I → R I 1 ∪I 2 such that we can write the above discrete compatibility condition as In order to approximate the continuous energy (6) with our discrete Raviart-Thomas field σ RT we can use different types of quadrature rules. After various attempts, we found out that simply summing the values at the center points of the volumes V i,j,k provides the highest flexibility for the discrete field σ to concentrate on thin lines, and, in turn, yields the most faithful numerical results.
Letting thus σ i,j,k = σ RT (iδ x , jδ x , kδ θ ), we have that this discrete field is obtained by the following formula: ). (24) We again introduce an operator A : R J 1 ∪J 2 ∪J θ → (R 3 ) J such that the above averaging operation can be writtenσ = Aσ.
Now, using the volume centered quadrature rule, the discrete energy is The consistency of this energy, as δ := (δ x , δ θ ) → 0, with the continuous energy F defined in (6) is studied in Appendix A.
The regularization function h(θ, p), where θ ∈ S 1 and p = (p x , p θ ) ∈ R 3 is defined as in (3), and the function f (t), which appears in its definition will be one of the following classical examples of convex functions: In all examples, the regularizing function is a combination of length and curvature regularization and the parameter α > 0 can be used to adjust the influence of the curvature regularization.
The first function, f 1 is the sum of length and absolute curvature, hence we call the corresponding regularizer "total absolute curvature" (TAC). One of its main features is that it allows for sharp corners in the level sets of the image. Interestingly, since integrating the absolute curvature along the boundary of a shape is constantly 2π for all convex shapes (and in general, is scale independent) then the main effect of the curvature term in this energy is to penalize non-convex shapes, regardless of their size.
The function f 2 combines length and curvature through an Euclidean metric and hence, it corresponds to the total variation of the lifted curve in the RT space, hence we consequently denote this regularizer "total roto-translational variation" (TRV). For relatively small curvature, it favors smooth shapes, but it also allows sharp discontinuities.
Function f 3 penalizes squared curvature plus length and hence is equivalent to the Elastica energy. We call this regularizer "total squared curvature" (TSC). It is well-known that while length regularization favors smaller shapes, quadratic curvature favors larger shapes. Hence, the interplay between length and curvature regularization removes the shrinkage bias and leads to smooth shapes.
Based on the "basis" functions f l , l = 1, 2, 3, we obtain the following convex regularization functions in the RT space.
Observe in particular, that the above definitions already properly take into account the correct values of the function h l , l = 1, 2, 3 for p x = 0. In the following sections we will apply the proposed curvature based regularization functions to a variety of image-and shape processing problems. For this, we consider generic optimization problems of the form min where G(u) is a convex, lsc. function defining an image-based convex data fidelity term, possibly dependent on the pixel location.

Primal-dual optimization
In this section, we show how to compute a minimizer of the non-smooth convex optimization problem (29). For notational simplicity, we first divide the whole objective function (29) by δ 2 x δ θ and assume that the remaining factors are absorbed by the data fitting term.
In order to solve the constrained optimization problem (29), we consider its Lagrangian (saddlepoint) formulation: where φ ∈ R J , ξ = (ξ 1 , ξ 2 , ξ θ ) ∈ (R 3 ) J and ψ = (ψ 1 , ψ 2 ) ∈ R I 1 ∪I 2 with ψ 1 ∈ R I 1 and ψ 2 ∈ R I 2 are the dual variables (Lagrange multipliers or discrete test functions). The function h * denotes the convex conjugate of the function h. Recall that h is the support function of the convex set where f * denotes the convex conjugate of f . Hence, the convex conjugate h * is simply the indicator function of the set H: The problem (30) is a saddle-point problem which is separable in the non-linear terms and hence falls into the class of problems that can be solved by the first-order-primal-dual algorithm with diagonal preconditioning and overrelaxation [24,60,25]. In order to make the algorithm implementable, we need efficient algorithms to compute the projection operators proj H(θ) and proximity operators prox τ G . They are defined as the unique minimizers of the minimization problems: proj H(θ) (η) = arg min and, for some τ > 0. In what follows, we will detail the projection operator for different instances of the convex sets H (respectively regularization functions h), the proximity maps for G will be detailed as soon as they are needed in the numerical results. The general idea for performing the projection ξ = proj H(θ) (η) of a point η = (η x , η) onto the set H(θ) is outlined in Figure 3. The Figure on the left hand side represents the set H in the ξ x plane, and the Figure on the right hand side shows the "profile" of the set which is obtained by cutting along the ξ x θ × ξ θ plane. We let η x θ = η x · θ and define (ξ x θ , ξ θ ) = proj P (η x θ , η θ ) as the projection of the point (η x θ , η θ ) onto the "profile" Then, using simple geometric reasoning, we see that the variable ξ x can be recovered from ξ θ by It remains to detail the projections onto different profiles P . Figure 3: Projection of the point (η x , η θ ) onto the convex sets H. First, the point (η x θ , η θ ) is projected onto the profile P which yields (ξ x θ , ξ θ ), then ξ x is computed from ξ x θ via the geometric relationship ξ • TAC: f 1 (t) = 1 + α|t|: The convex conjugate of the function f 1 is given by and in turn the profile P 1 is given by It is straightforward that the projection of a point (η x θ , η θ ) onto the profile P 1 can be performed via simple truncation operations .
• TRV: f 2 (t) = √ 1 + α 2 t 2 : A simple computation shows that the convex conjugate of f 2 (t) is given Inserting the expression of the convex conjugate into the profile (31) and squaring both sides, we obtain In what follows, we assume that max(0, η x θ ) 2 + (η θ /α) 2 > 1 since otherwise we do not need to project the point. We first treat the case η x θ ≤ 0. It is easy to see that in this case the solution of the projected point is given by min(α, η θ )) .
In case η x θ > 0, computing the projection of a point (η x θ , η θ ) onto the boundary of P 2 amounts to solve the following equality constrained optimization problem The Karush-Kuhn-Tucker (KKT) optimality conditions for the above problem are given by where λ > 0 is the Lagrange multiplier which is positive since the point (η x θ , η θ ) was assumed to be outside P 2 . Combining the first three equations shows that the Lagrange multiplier is computed from the roots of the fourth-order polynomial Let us observe that in case α = 1, the solution for λ is particularly simple, indeed such that from the first two equations we obtain for the projected point In the general case α = 1 we compute λ by applying Newton's algorithm to find the (correct) root from the fourth order polynomial. In our experiments it turns out that a large enough initial value, e.g. λ = 10 3 provides a good initialization for Newton's algorithm. Usually, we need less than 5-10 iterations of Newton's algorithm to converge to a solution with feasibility error less than 10 −9 . From the first two equations of the KKT conditions the expression of the projected point is given by: .
• TSC: f 3 (t) = 1 + α 2 t 2 : The convex conjugate of f 3 is computed as and hence the profile P 3 is given by Following the same approach as before a point (η x θ , η θ ) with η x θ + (η θ /(2α)) 2 > 1 is projected onto P 3 by solving the KKT optimality conditions Combining the above equations, the Lagrange multiplier is obtained by finding the correct root from the third-order polynomial We again apply Newton's method and observe that for a large enough initial value, e.g. λ = 10 3 , Newton's algorithm rapidly converges to the correct root of the polynomial. The projected point is finally given by

Computing a disk
In the first example, we consider the most basic numerical experiment, which is using TSC energy to compute a disk of a given radius. The aim of this experiment is to investigate the quality of our proposed discretization scheme, in particular when using a different number of discrete orientations. Consider a 2D disk B(0, r) = {x ∈ R 2 : |x| ≤ r} of radius r > 0 centered around the origin. The TSC energy of the boundary of the disk B(0, r) is given by Here our disk will be represented by its characteristic function χ B(0,r) (x). We can force the minimizer of the TSC energy to yield a disk using suitable boundary conditions. For example, we can force at least one point inside the disk to be one and at least one point outside the disk to be zero. The following simple computation shows that the minimizer of the TSC energy will be a disk of radius r = α, indeed: which is minimized if κ = 1/r = 1/α everywhere. We set up a our optimization problem using a 2D grid of N 1 × N 2 = 40 × 40 pixels. We will use a different number of discrete orientations N θ ∈ {4, 8, 16, 32, 64} in order to investigate the quality of the approximation depending on N θ . In this and the subsequent experiments, the discretization width of the spatial grid is set to δ x = 1 and the discretization width of the angular dimension is set to δ θ = 2π/N θ . The aim of this experiment is to compute a disk of radius r = 10 pixels. We therefore define an inpainting domain D ⊆ I forming a band of 10 pixels width around a disk of radius r = 10. In this domain we minimize the TSC energy and we use the remaining part of the image as boundary condition. The corresponding data term in our optimization problem (29) is given by where ι C denotes the indicator function of the convex set C. The proximal map for G(u) is given by  In order to quantify the approximation quality of our discretization for a varying number of discrete orientations, we report the total variation (TV), the absolute curvature (AC) and the squared curvature (SC). All three quantities are computed from the averaged fieldσ = Aσ using where the function h is one of the three instances: . Table 1 details the values of the discrete energies we obtained for different numbers of discrete orientations N θ . From the results, one can see that a higher number of discrete energies generally leads to a better approximation of the true energy. This is particularly true for the value of the squared curvature which seems to be well approximated only when using a quite high number of discrete orientations. The absolute curvature, however seems to be well approximated even when using only a small number of discrete orientations. The reason for this is that the absolute curvature of a smooth curve is easy to approximate by means of a piecewise linear curve. Figure 4(f) visualizes the averaged Raviart-Thomas vector fieldσ of the disk example in case of N θ = 64 discrete orientations. Observe that the vector field nicely corresponds to the expected shape of a helix, shown in Figure 1. However, due to diffusive effects of our numerical scheme, the vector field does not perfectly concentrate on a one-dimensional structure.

Non-smooth level sets
Here we demonstrate a typical effect of our convexification which finds a low-energy solution for an image with non-smooth level sets which should have infinite energy in more standard relaxations [53] of the Elastica energy. Figure 5(a) shows an image of size N 1 × N 2 = 100 × 50 pixels of a black and a white disk in front of a gray background. Similarly to the inpainting problem of the previous example, we fix only the four small parts of the original image (see Figure 5(b)) and minimize the TSC energy in the inpainting domain D, indicated by the gray area. In this experiment we used N θ = 64 discrete orientations and we set α = 17 to match the radius of the disks. Figure 5(c) shows the computed minimizer of the TSC energy. Observe that the solution u does not yield the expected two disks of the original image but rather drop-like shapes that form a sharp cusp in the middle of the image. Figure 5(b) shows a stream-line representation of the  minimizing vector fieldσ in the roto-translation space. Inspecting the field, one can immediately see the reason for this behavior. The vector fieldσ forms a twisted ∞-shape curve which "skips" the strong curvature of the cusp.

Shape completion
In this section we provide some qualitative results on a number of different shape completion problems. Similar to the previous two examples, we define an inpainting domain D which is indicated by the gray area and keep the remaining image as boundary condition. Figure 6 shows various input shapes with their inpainting domains and the solutions of minimizing different curvature energies using different settings of the parameter α. One can see that while minimizing the TAC energy usually leads to straight connections and sharp corners, TSC leads to a smooth continuation of the boundaries. The TRV energy leads to results which are somewhere in between the results of TAC and TSC. In the last row of Figure 6, we additionally demonstrate the behavior for completing straight lines at different rotations by minimizing the TAC energy. It turns out that  : Shape completion for various shapes using total absolute curvature (TAC), total rototranslational variation (TRV) and total squared curvature (TSC). The inpainting domain is indicated by the gray pixels. In the last row we provide results for completing straight lines at different angles.
our discretization scheme performs quite well for rotations of 0, π/4, π/2, . . . but leads to more diffusive results for rotations of π/8, 3π/8, . . .. The development of a more isotropic scheme will be subject of future research. In Figure 7, we show the application to shape completion from dipoles using the original data of "Weickert's cat" [76]. The input image u 0 is of size N 1 × N 2 = 400 × 380 pixels and the inpainting domain D is again indicated by the gray pixels. In this example we used N θ = 64 discrete orientations and the curvature parameter was set to α = 50. In order to avoid relaxation artifacts (self-intersections) inside the dipoles, we additionally use zero-boundary conditions for the field σ 1,2 within the constant areas of the dipoles. The results show that our numerical scheme can successfully reconstruct the curvilinear shape of the cat from only very little information given by the dipoles. However, we can also observe diffusion artifacts of our numerical scheme, especially if the reconstructed boundaries are relatively long. Sharper results for this kind of problems are usually obtained using sophisticated anisotropic diffusion schemes such as the edge enhancing diffusion (EED) [75,67]. It would be interesting to understand whether these schemes also minimize an underlying variational energy.

Shape regularization
In our next experiment, we apply our curvature based energies for shape regularization. Given an input image (u 0 i ) i∈I , which is the characteristic function of a given shape, our aim is to compute a simplified (or regularized) shape which is represented by means of a (relaxed) binary image (u i ) i∈I . We make use of a simple linear fidelity term which is frequently used in image segmentation: where w ∈ R I is a force field. For the application to shape regularization we use w = λ( 1 2 − u 0 ), where λ > 0 defines the strength of the force field. The proximal map for this data term is easily computed: In Figure 8, we apply TAC, TRV and TSC regularization to regularize the shape of Picasso's "Bull fight" image. In all three cases, we set α = 10 and we use different settings of the parameter λ to obtain gradually simplified shapes. From the results one can see that TAC regularization yields shapes with relatively straight boundaries and sharp corners. TRV regularization yields smooth shapes but also allows for sharp corners. TSC regularization yields smooth shapes. Observe that whenever it seems energetically preferable, the solution of TSC produces sharp cusps with "hidden" edges to bypass locations of strong curvature, as predicted by the theory [12]. This effect is usually less visible when minimizing the TAC or TRV energies.

Image inpainting
In this section, we apply our proposed curvature energies to the classical problem of image inpainting. Similar to the previous examples, we use the data term (32) with the only difference that the input image u 0 is now a gray level image.
In both examples we used N θ = 32 discrete orientations and the curvature parameter was set to α = 15. Figure 9 shows the results for two different inpainting problems. In the "Einstein" image, we randomly remove 80% lines and in the "Picasso" image we randomly remove 90% of the pixels. In case of the "Einstein" image, we use TSC regularization with α = 10 and for the "Picasso" we set α = 5. In case of the "Einstein" image, the gaps a much larger and hence we use a larger value of α. For comparison, we also provide results of standard TV regularization, which is equivalent to using α = 0 in one of the three curvature energies. From the results one can clearly see that curvature regularization leads to significantly better inpainting results. On the downside, we can also observe some artifacts which are caused by our convex representation in the roto-translation space.

Image denoising
In our last experiments we investigate our curvature energies for the classical problem of image denoising. We investigate two different types of noise: Zero-mean Gaussian noise and impulse noise such as "salt & pepper" noise. For Gaussian noise, it is well-known data a quadratic 2 data term where λ > 0 is a good data fidelity parameter. The proximal map is given by In case of impulse noise, a 1 data term is more suitable: since it is more robust with respect to outliers. The proximal map for the 1 data term is given by the classical soft-shrinkage formula In Figure 10 we show the results of our curvature-based regularization energies for image denoising. In case of the "Louvre" image we generated the noisy image by adding 25% "salt & pepper" noise. The noisy "Leberblümchen" image was generated by adding zero-mean Gaussian noise with standard deviation 0.1. In case of Gaussian noise we use a 2 data term and in case of "salt & pepper" noise we used the 1 fidelity term. In both cases, we use TSC regularization and the curvature parameter was set to α = 10. For comparison, we also show the results of standard total variation (TV) denoising, which is obtained from our models by setting α = 0. One can clearly see that TSC regularization leads to a better preservation of image edges than TV regularization, in particular at small and elongated structures. This effect is more visible in case of "salt & pepper" noise since the 1 data term leads to a behavior similar to image inpainting, once an outlier is detected.

Proof of Theorem 1
This section and the following are devoted to the proof of Theorem 1. We first show a preliminary result which shows that admissible curves (of finite energy) in Ω × S 1 project onto curves in Ω with bounded energy as well. and assume that it is finite, which in particular yields that for a.e. s, γ x (s) = |γ x (s)|(cos γ θ (s), sin γ θ (s)) T .

Control of curves
We consider the projected curve γ x : [0, L] → Ω (possibly overlapping even if γ is simple) and wish to show that the energy E controls an energy on this curve. First we reparameterize (as usual) the rectifiable curve γ x by its length. For this we introduce the length, for s ∈ [0, L] ( (s))θ( (s))˙ (s)ds which shows that θ(l) is a BV function, or equivalently thatγ has a curvature κ =θ which is a bounded measure; in particular taking the supremum over ψ with |ψ| ≤ 1 one finds In addition, it follows, given ψ ∈ C ∞ c ( (a), (b)) with 0 ≤ a < b ≤ L that: and taking the supremum one deduces that We now show the following lemma: Lemma 5.1. Let Γ ⊂ Ω be a C 2 (oriented) curve with tangent τ Γ and curvature κ Γ , ane let and γ : [0, L] → Ω × S 1 be a rectifiable curve. Define Γ + = {x ∈ Γ : ∃s ∈ (0, L), γ x (s) = x andγ x (s) · τ Γ ≥ 0}. Then for any bounded, nonnegative, continuous function ϕ ∈ C 0 (Ω; R + ), Proof. We just need to show that, defining the measure κ as before, and the conclusion will follow from (34).
A first observation is that as τ (l) defined above is BV , it has at most a countable number of jumps and one can cover (up to the jump points) [0, (L)] with an at most countable union of intervals (a i , b i ) i∈I on which τ (l) is continuous (henceγ is a C 1 curve on (a i , b i )). Moreover, possibly dropping further a finite number of points, one may assume that |κ|(a i , b i ) < π so that γ x (a i , b i ) is a simple curve.
Assume |{l ∈ (a i , b i ) :γ(l) ∈ Γ,γ(l) · τ Γ ≥ 0}| > 0. Choose such an l ∈ (a i , b i ), with x =γ(l) ∈ Γ + , and such thatγ is differentiable in l. We assume in addition that x is a point of (H 1 -)density one inγ(a i , b i ) ∩ Γ + , that l is a Lebesgue point of κ a , the absolutely continuous part (w.r. the Lebesgue measure) of κ, and that θ(l) is approximately differentiable at l: for any η > 0, All this is true H 1 -a.e. inγ(a i , b i ) ∩ Γ + , so there is no loss of generality, see for instance [3]. Ifγ(l) = τ Γ (x), then x must be an isolated point inγ(a i , b i ) ∩ Γ + , a contradiction. Hencė γ(l) = τ Γ (x). In the same way, one has that for a.e. s near l such thatγ(s) ∈ Γ,γ(s) = τ Γ (γ(s)). For such s, one therefore has that where θ Γ is the angle between 1 0 and τΓ. Using that for such s, we obtain that for such s, (where o(1) goes to 0 as ε → 0). Hence, for ε > 0 small and η > 0, In the limit, using that l is a point of density one in the set {γ(·) ∈ Γ}, we obtain that |κ Γ (x) − κ a (l)| ≤ η and since η is arbitrary, κ Γ (x) = κ a (l). Since the intervals (a i , b i ) cover H 1 -almost all of [0, (L)], we find that this equality holds H 1 -a.e. in Γ + . The thesis of the Lemma easily follows, since

Decomposition of σ
Let u ∈ BV (Ω) with F (u) < ∞. We observe that we can decompose an admissible σ (with in particular div σ = 0) as follows: where µ is a nonnegative measure on the set of charges M which is here a shorthand notation for M(Ω × S 1 ; R 3 ): see Appendix B which discusses the results of Smirnov in [73].
Corollary 5.3. Let σ be admissible for problem (6), and assume it is decomposed as in (36).
Corollary 5.4. Let σ be admissible for problem (6), and assume it is decomposed as in (36). Let h be defined as in (3), and ϕ ∈ C 0 c (Ω) a nonnegative test function. Then Proof. This is a consequence of Lemma 5.2, however now h can take the value +∞. We simply observe that if we let, for ε > 0, then h ε is (1/ε)-Lipschitz, convex, one-homogeneous in the second variable, and continuous (and hence also (ξ, p) = ((x, θ), p) → ϕ(x)h ε (θ, p)). Hence (37) yields The results follows from the monotone convergence theorem.

Conclusion: proof of Theorem 1
We can now complete the proof of Theorem 1. We consider E a C 2 set. First, if γ(t) : [0, 1] → Ω is an oriented parameterization of (the closure of) ∂E ∩ Ω, (assuming it is connected, otherwise one needs to introduce a curve for each connected component), one can define in Ω × S 1 the curve λ(t) = (γ(t), θ(t)) where θ(t) is defined by λ(t) = |γ (t)|(cos θ(t), sin θ(t)) T . We then let Γ = λ([0, 1]) ∩ Ω × S 1 and σ = τ Γ H 1 Γ. Then, σ is admissible for (6) with u = χ E , and one finds that We need therefore to prove the reverse inequality. Let σ be admissible for (6). Thanks to [73, Theorem A] (cf Appendix B, eq. (59)), one can decompose σ as where λ are of the form λ γ = τ γ H 1 γ for rectifiable (possibly closed) curves γ ⊂ Ω × S 1 of length at most one. Equation (39) is valid for this decomposition and shows that for any ϕ ∈ C 0 c (Ω; [0, 1]), Now, for any ψ ∈ C 0 c (Ω; R 2 ), one has where τ E is a tangent vector to ∂E (oriented with E on the left-hand side and E c on the right-hand side 1 ). Let us introduce the signed distance function d E (x) = dist (x, Ω \ E) − dist (x, E) which is C 2 in a neighborhood of ∂E and is such that ∇d E = ν E (the inner normal) on ∂Ω and ∆d E = −κ E (we assume the curvature is nonnegative where the set is convex). If we consider, for ε > 0, a test function of the form Sending ε → 0 we find that Now, we can choose ϕ of the form ϕ(x)f (κ E (x)), with ϕ ∈ C 0 c (Ω; [0, 1]), assuming f is finite-valued. We obtain that Given the rectifiable curve γ, one has that thanks to (35) in Lemma 5.1. Then, (40), (42) and (43) yield and Theorem 1 is easily deduced. In case f takes the value +∞, then we first approximates f with a finite-valued function from below (replacing f with min t f (t ) + |t − t |/ε for ε > 0 small), and once the inequality is established for this function we send ε → 0, so that (9) also holds.

A Consistency of the discretization
In this appendix, we study the consistency of the discrete approximation of the problem which is used in Section 4.

A.1 Preliminary results
Let us introduce, for (s, t) ∈ R 2 ,h which is such thath (s, t) = sup as + bt.
Proof. If we introduce the (averaged) marginal σ ∈ M 1 (S 1 ; θ) defined by S 1 ψσ = 1 δ 2 x S×S 1 ψσ for all ψ ∈ C 0 (S 1 ), then one observes that This follows by a disintegration argument and using Jensen's inequality in each "slice" corresponding to a fixed value of θ. The result then follows from (44), together with Remark A.1. The last statement comes from the fact that σ x is the average of measures all contained in the cone, which is convex.
Consider now a measure σ admissible for some function u, and assume that σ is "smooth" in x: we assume for instance that it is the result of a convolution ρ ε * σ for some σ admissible (possibly extended in a larger domain), with ρ ε (x) a rotationally symmetric mollifier, as in the proof of Proposition 3.1. In this case, x → σ can be seen as a M 1 (S 1 ; R 3 )-valued smooth function.
Consider (x,θ) ∈ Ω × S 1 , δ x , δ θ small, and define in the volume V = S × [θ − δ θ /2,θ + δ θ /2) the averageσ as before and the "Raviart-Thomas" approximation of σ defined by the average fluxes through the 6 facets of V (linearly extended inside the volume, as in eq. 21). There are several ways to define this properly, at least for allθ but a countable number. In our case, one can disintegrate the measure in Ω × S 1 as σ = σ θ dµ where µ is a bounded positive measure in S 1 and for all θ, σ θ ∈ C ∞ (Ω). Then From Corollary A.5, letting e i,j,k =σ i,j,k −σ i,j,k , one has δ 2 x δ θ Proposition A.7. Let u ∈ BV (Ω), σ be admissible for u and such that Then for δ = (δ x , δ θ ) → 0 one can find a discrete field (σ 1 ) and a discrete image u δ (strongly in L 2 (Ω)) and such that for all i, j, (23) holds, for all i, j, k, (51) and (22) hold, and: whereσ is defined by (24), and whereh δ is defined in (50) (or ish in case it is Lipschitz).
To show that the discretization is consistent, we must now show a similar lower bound: namely that given any u and u δ , σ,σ which satisfy (24), (51), (23), and (52) (weakly, for instance as distributions), then one has A first obvious remark is that the field is bounded in measure, and hence, up to subsequences, converges (weakly- * ) to a measure σ. It is then easy to deduce from (51) and the convexity ofh that One can also check that div σ = 0 by passing to the limit in (22) (after a suitable integration against a smooth test function, exactly as in (55) below). Hence it is enough to show that the limiting σ is compatible with u. But this is quite obvious from (23), which one can integrate against a smooth test function, then "integrate by part" before passing to the limit. More precisely, for ϕ ∈ C ∞ c (Ω), one has (dropping the superscripts δ and denoting ϕ i,j = (1/δ 2 x ) Si,j ϕ(x)dx):  (55) as δ → 0. Eventually, we need to show a compactness property, which is that if sup δ δ 2 x δ θ i,j,kh (σ i,j,k · θ k ,σ θ i,j,k ) < +∞ (56) the discrete image u δ which is recovered from (23) (up to a constant) converges to a u(x), x ∈ Ω (in a weak sense which will be made clear). The point here is that a priori, from (56) and (1) < +∞ (57) so that the discrete total variation of u δ is a priori not well controlled. However, one can easily check that the kernel of the operator which appears in the energy (57) is two-dimensional, and made of the oscillating discrete images α, β ∈ R 2 . Hence it is possible to show that one can decompose u δ as a sum of a non-oscillating function with zero averageū δ and an oscillation v δ , and obtain a strong control on the discrete total variation ofū δ . Therefore one easily deduce that any suitably built continuous extension of u δ will converge to some u strongly in L p (Ω), for any p < 2 (as BV (Ω) is compactly embedded in such spaces), and weakly in L 2 (Ω). In addition, any control on the average of u δ and on its oscillation (which cannot be given by (57) and has to come from other terms in the energy, such as a boundary condition or a penalization: note that it is enough to control two adjacent pixels) will ensure in addition that v δ remains bounded and converges (only weakly in L p (Ω), if the control is only on the L p norm, however in this case it is obvious that the oscillating term in (58) goes to zero and v δ can only go to a constant).
To sum up, we have shown the following.
Remark A.9. In practice, we did not use the inf-convolutionsh δ (onlyh) in our discrete scheme. Also, we replaced the constraint (51) with the stronger constraintσ x i,j,k ∈ R + θ k , after having experimentally observed that there was no qualitative difference in the output. It seems the results we compute are still consistent with what is expected from the energy. B Smirnov's theorem in Ω × S 1 In this whole paper Ω ⊂ R 2 is assumed to be a Lipschitz set. In particular, locally its boundary can be represented as the subgraph {(x, y) : y < h(x)} of a Lipschitz function h. Consider a ball B where this representation holds and assume first h is C 1 , then one can extend in B a bounded Radon measure σ with div σ = 0 intoσ defined (for ψ ∈ C 0 c (B; R 2 )) Bσ · ψ := B∩Ω σ · ψ(x, y) − 1 2h (x) 0 −1 ψ(x, 2h(x) − y) .
A similar construction would allow to extend a field σ ∈ M(Ω × R 2 ; R 4 ) to M(R 4 ; R 4 ) with free divergence (either in a neighborhood or Ω × R 4 , or even everywhere). This remark allows to localize Smirnov's theorems in [73].
Consider indeed now a free divergence field σ ∈ M(Ω × S 1 ; R 3 ). It can be seen, after extension, as a field in M(R 2 × R 2 ; R 4 ) with sptσ ⊆ R 2 × S 1 .
As in Smirnov's paper [73], for l > 0 we introduce C l the set of oriented curves γ in R 4 with length l, with the topology corresponding to the weak convergence of the measures τ γ H 1 γ. Then, thanks to [73,Theorem A], σ can be decomposed as σ = C l λdµ(λ), |σ| = C l |λ|dµ(λ), for some measure µ on C l . Moreover thanks to Remark 5 in [73], µ-a.e. curve in the decomposition lies in R 2 × S 1 .
For this work, it is enough to consider l = 1. Moreover, if we restrict then all these measures to Ω × S 1 (and take for µ the corresponding marginal), we get a decomposition on curves of length less or equal to 1 (possibly entering/exiting the domain). By a slight abuse of notation we still denote C 1 such a set of curves. One finds that with now λ ∈ C 1 , the curves of length at most one in Ω × S 1 .
Remark B.1. Theorem B in [73] is a more precise statement. It shows that one can obtain a similar decomposition with now curves λ with div λ = 0 a.e.: being either finite curves entering and exiting the domain, or "elementary solenoids", which are objects of the form (in particular one should have |f (t)| = 1 a.e.), meaning that for any ϕ ∈ C 1 c (Ω × S 1 ), This expresses that either λ is defined by the closed curve f (R) (if f is periodic), or λ is a limit of curves which densify and do not loose mass in the limit. We do not need such a precise result for our construction.