ITERATIVE FINITE ELEMENT SOLUTION OF A CONSTRAINED TOTAL VARIATION REGULARIZED MODEL PROBLEM

. The discretization of a bilaterally constrained total variation minimization problem with conforming low order ﬁnite elements is analyzed and three iterative schemes are proposed which diﬀer in the treatment of the non-diﬀerentiable terms. Unconditional stability and convergence of the algorithms is addressed, an application to piecewise constant image segmentation is presented and numerical experiments are shown.

1. Introduction. In this paper we investigate iterative schemes for the solution of the bilaterally constrained ROF problem which consists in the minimization of the functional E(u) := |Du|(Ω) + α 2 u − g 2 L 2 (Ω) + I K (u) in the set of functions with bounded total variation |Du|(Ω), where I K is the indicator functional of a convex set K ⊂ L 2 (Ω) characterized by two obstacle functions. While existence and uniqueness of a minimizer of E can be established using the direct method in the calculus of variations, the non-differentiability of the total variation and the indicator functional are challenging from a numerical point of view. The minimization problem serves as a model problem for a wide class of functionals involving the total variation and the indicator functional which comprises a bilateral constraint on the variable u. For example, in [6,14] the above minimization problem is considered in the context of image denoising and the bilateral constraint is imposed in order for the restored image to have pixel values lying in a certain range, e.g., in the interval [0,255]. In [13] Chan et al. show the equivalence of a minimization problem occuring in image segmentation to the minimization of a functional that resembles E with the quadratic fidelity term replaced by a linear functional in u. The application of a regularization term to a bilaterally constrained variable also occurs in the modelling of damage processes in continuum mechanics, where the bilaterally constrained variable is the damage variable and the regularization term consists of a differentiable L p -norm (p ≥ 2) of the gradient of the damage variable (cf. [19]) or the total variation of the damage variable (cf. [21]) while the constraint ensures that the damage variable takes on values only in the interval [0, 1].
In the recent years many PDE-based algorithms have been proposed for the classical ROF functional proposed by Rudin, Osher and Fatemi in [20], where the bilateral constraint is omitted, reaching from primal-dual methods, combinations of regularization and standard minimization techniques and dual methods to the alternate direction method of multipliers, augmented Lagrangian methods and split Brègman methods, see, e.g., [3,23] for an overview of many of those algorithms. Regarding the bilaterally constrained ROF problem Casas et al. proposed in [10] an active set strategy using the anisotropic BV -seminorm. This ansatz works with the primal formulation of the problem and does not use regularization. However, a penalization approach is used to solve auxiliary subproblems. In [6], motivated by the dual approach of Chambolle in [11], Beck and Teboulle suggest a monotone fast iterative shrinkage/thresholding algorithm (MFISTA) for the dual of the bilaterally constrained ROF problem for which the convergence rate O(1/k 2 ) with k being the iteration counter is proven. Yet, as in [11], the constant entering the convergence rate grows like O(1/h 2 ) with h being the mesh size of the underlying grid. Chan et al. proposed in [14] an augmented Lagrangian method, which has already been successfully applied to the ROF problem by Wu and Tai in [22]. Let us remark that the primal-dual algorithm proposed in [12] can also be applied to solve the bilaterally constrained ROF problem.
The aim of this paper is to compare three different approaches for the iterative minimization of E and to identify the most appropriate method. To this extent, we introduce a finite element discretization with globally continuous piecewise affine finite elements. The iterative algorithms differ in the treatment of the two challenging terms appearing in E. The first approach uses splitting techniques for both the total variation and the bilateral constraint and is a variant of the augmented Lagrangian method proposed in [14], where the involved norms are properly weighted in order to ensure unconditional stability of the method, see also [3] where the authors proposed a properly weighted augmented Lagrangian method for the classical ROF problem. The next ansatz uses regularization to handle the non-differentiability of the total variation and the indicator functional and employs a semi-implicit subgradient flow which turns out to be unconditionally stable and will be referred to as the Heron-penalty method. For the classical ROF problem the Heron method has already been introduced in [3]. The third algorithm uses a combination of the regularization and the splitting approach, i.e., using the observation made in [3] that a regularization approach seems to be the most appropriate to deal with the total variation, we employ a regularization for the total variation and a splitting technique for the bilateral constraint, which we call the Heron-split method. The primal-dual algorithm proposed in [12] has been compared to the Heron method and the augmented Lagrangian method in [3] for the classical ROF problem and it has been observed that the primal-dual algorithm is slow due to the restrictive condition on the step size.
The outline of the paper is as follows. In Section 2 we give some basic notation and introduce the space of functions with bounded variation as well as the globally continuous piecewise affine and the piecewise constant finite element space. We present the minimization problem in Section 3 and discuss existence and uniqueness of a minimizer. In Section 4 we consider the finite element discretization and prove convergence of discrete minimizers to the continuous solution before we present the three algorithms and address their stability in Section 5. In Section 6 we address the minimization of the functional prove convergence of discrete energies to the minimal energy and recall a minimization problem arising in an image segmentation model discussed by Chan et al. in [13]. We finally provide results of numerical experiments for the three iterative schemes applied to the bilaterally constrained ROF problem and for a particular choice of g and two different values of α as well as for the image segmentation problem in Section 7. We then conclude the paper in Section 8.
For more details concerning the space BV (Ω) we refer, e.g., to [1]. For a convex functional F : X → R ∪ {∞} on some Hilbert space X with inner product (·, ·) X we define its subdifferential ∂F (u) at u ∈ X by If F is Fréchet-differentiable at u, we denote by δF (u) its Fréchet-derivative at u. For a sequence (a j ) j∈N and a step size τ > 0 we let d t a j+1 := a j+1 − a j τ be the backward difference quotient.

2.2.
Finite element spaces. Throughout the paper we let c be a generic constant. We denote by (T h ) h>0 a family of uniform triangulations of Ω with mesh sizes h = max T ∈T h h T where h T is the diameter of the simplex T . For a given triangulation T h the set N h contains the corresponding nodes and we consider the finite element spaces of continuous, piecewise affine functions and of elementwise constant functions (r = 1) or vector fields (r = d) We also denote by T , ∈ N, a triangulation of Ω generated from an initial triangulation T 0 by uniform refinements and the integer being related to the mesh size h by h = c2 − . The set of nodes N is defined correspondigly.
We introduce on S 1 (T h ) the discrete norm · h induced by the discrete scalar product is the nodal basis function associated with the node z ∈ N h . The norm · h and the L 2 -norm are related on S 1 (T h ) as follows.
3. Model problem. For a function g ∈ L ∞ (Ω) and a parameter α > 0 we consider the constrained minimization problem defined by the functional in the set of functions BV (Ω) ∩ L 2 (Ω). The convex set K is given by with χ, ψ ∈ BV (Ω) ∩ L ∞ (Ω) and χ < ψ almost everywhere. We will refer to this minimization problem as the bilaterally constrained ROF problem.
Proof. The set of feasible functions is nonempty and the functional E is bounded from below. Hence, there exists an infimizing sequence (u j ) j which is bounded in BV (Ω) ∩ L 2 (Ω) due to the coercivity of E. The compact embedding of BV (Ω) into L 1 (Ω) (cf. [1]) and the boundedness of the sequence imply the existence of anot relabelled -weakly convergent subsequence u j u in BV (Ω) and by passing to a further subsequence u j u in BV (Ω)∩L 2 (Ω). Since u j → u in L 1 (Ω) there exists a further subsequence such that u j → u almost everywhere in Ω. This implies that χ ≤ u ≤ ψ almost everywhere in Ω since χ ≤ u j ≤ ψ for all j ∈ N. The weak lower semicontinuity of E (cf. [1]) then yields that u is a minimizer of E. Uniqueness follows by the strict convexity of E on K.
We define the functionals G : Lemma 3.1. Let u be the unique minimizer of E in BV (Ω) ∩ L 2 (Ω). Then for any v ∈ BV (Ω) ∩ L 2 (Ω) we have Proof. We note that G is Fréchet-differentiable, i.e., we have Standard arguments from convex analysis (cf. [1, Thm. 9.5.4]) yield The strong convexity of G and the subgradient inequality then imply that for arbi- which proves the lemma.

4.
Finite element discretization. We will need the following assumption for the error analysis which imposes restrictions on admissible domains and obstacles.
We seek a discrete minimizer u h ∈ S 1 (T h ) of where we used that |Dv|(Ω) = ∇v L 1 (Ω) for v ∈ W 1,1 (Ω). Existence and uniqueness can be proven following the arguments of the proof of Proposition 1.

Proposition 2.
Let u and u h be the unique minimizers of E in BV (Ω) ∩ L 2 (Ω) and S 1 (T h ), respectively. Then we have for d = 2 Proof. Let d = 2. Following the arguments in the proof of [2, Thm. 10.7] we note that by Lemma 3.1 we have for v = u h and for arbitrary v h ∈ S 1 (T h ) that Choose The choice ε = h 1/2 yields the first estimate. If d = 3, we proceed as above and obtain by Lemma 4.1 The choice ε = h 1/3 gives the second estimate.
Remark 3. The error estimate is suboptimal in the sense that for u ∈ BV (Ω) ∩ L ∞ (Ω) we have the approximation property 5. Iterative solution. The iterative minimization of E is difficult due to the nondifferentiability of the seminorm and the occurrence of the constraints. In this section we will discuss three approaches to deal with these difficulties.
to enforce these relations. Note that since the sequence (∇u h ) h>0 of discrete minimizers u h of E may not be bounded in L 2 (Ω) but only in L 1 (Ω) we have to measure these quantities in a weaker norm than the L 2 -norm. We therefore equip the space and induced norm · w . An inverse estimate shows that (∇u h ) h>0 is bounded with respect to · w , cf. [3]. We then consider the consistently stabilized Lagrange functional The parameters σ 1 and σ 2 are assumed to be positive. We then have that (u h , p h , [22,Thm. 4.1] for the same result for the ROF functional without bilateral constraint). Indeed, since cf. [18,Thm. 23.9], where the operator div h : for all (µ h , ν h ). One can show that the function . Using (2) and the particular choice of functions λ h and η h we observe that the triple (u h , p h , s h ) satisfies these optimality conditions. Altogether, we have Due to the fact that p h is elementwise constant and that mass lumping is used for the Lagrange multiplier η h and the auxiliary variable s h the minimization with respect to p h and s h can be done explicitly on each element and in each node, respectively. We approximate a saddle-point with the following iterative scheme (cf. [15]).
and stop if 1 stop . Increase j → j + 1 and continue with (1) otherwise.
Note that step (3) in Algorithm 1 resembles ascent steps for L h with respect to λ h and η h . The following stability estimate for the split-split method is a consequence of general results for splitting methods, see, e.g., [15,16], and arguments exploiting only the (strong) convexity of the involved functionals, which have been used by Wu and Tai in [22] for the classical ROF problem. It implies the convergence For arbitrary initializations of Algorithm 1 we have for any J ≥ 0

Regularization and penalization.
Another way to approximate the discrete minimizer u h of E is to regularize the BV -seminorm, introduce an auxiliary variable s h = u h and to penalize the difference s h − u h with a penalization parameter δ > 0 while allowing u h to penetrate the obstacles. We obtain the functional where |∇v h | ε = (|∇v h | 2 + ε 2 ) 1/2 and ε, δ > 0. The introduction of s h decouples the two nonlinearities. As in the split-split method, we use mass lumping for the variable s h to update it node-wise. A similar approach has been used by Hintermüller et al. in [17] for the predual of the classical ROF problem and an alternate minimization technique has been employed for the numerical solution.
Remark 4. Note that in contrast to the other two methods, the constraint u h ∈ K is enforced only in the limit as δ → 0.
We can establish the existence of a discrete minimizer for v h ∈ S 1 (T h ) and all z ∈ N h . The optimality conditions for a minimizing pair (u h,ε,δ , s h,ε,δ ) then imply the relation s h,ε,δ = Π K,h u h,ε,δ .
If T h consists only of right triangles (d = 2) or of tetrahedra having three mutually orthogonal edges (d = 3) and if both χ and ψ are globally constant, then we have for any T ∈ T h Proof. Consider d = 2, the case d = 3 is analogous. Without loss of generality we assume that the two edges conv{z 0 , z 1 } and conv{z 0 , z 2 } of T forming the right angle are parallel to the x-axis and the y-axis, respectively. Otherwise, we rotate the triangle which preserves the length of the gradients. We Let χ ≡ 0 and ψ ≡ 1 and define v h to be affine-linear on be the unique minimizer of E and a minimizer of E h,ε,δ , respectively, and let and particularly Proof. The proof of (3) is analogous to the proof of Lemma 3.1. Since (3). Observing that |x| ≤ |x| ε ≤ |x| + ε for x ∈ R d , and using Lemma 2.1 and the minimality property of u h , we deduce Using the identity s h,ε,δ = Π K,h u h,ε,δ , the triangle inequality, an inverse estimate, the binomial formula a − c 2 − b − c 2 = (a − b)(a + b − 2c) and Young's inequality we deduce which implies the second estimate. To prove the third estimate we note that similarly to Lemma 3.1 we have where the third inequality is due to the fact that By the second estimate we have 1 4δ s h,ε,δ − u h,ε,δ 2 ≤ c(ε + δh −2 ) and this implies the third estimate. Remark 7. We can significantly relax the restriction on δ if both χ and ψ are globally constant and if T h consists only of right triangles or of tetrahedra having three mutually orthogonal edges. Indeed, by Lemma 5.1 we in this case have ∇Π K,h u h,ε,δ L 1 (Ω) ≤ ∇u h,ε,δ L 1 (Ω) so that we do not need to use an inverse estimate in the proof of Proposition 4. This yields the improved estimates In this case, we may choose δ = h 1/2 /α for d = 2 (for d = 3: δ = h 1/3 /α) and δ = h/α to recover the convergence rate from Proposition 2 and the optimal convergence rate, respectively.
To define a stable scheme for the approximation of a minimizer (u h,ε,δ , s h,ε,δ ) ∈ Note that for fixed (v h , r h ) it follows that with optimal q h given by q h = (|∇v h | 2 + ε 2 ) 1/4 .
Due to the separate convexity of the functional E h,ε,δ with respect to v h , q h and r h we use a decoupled semi-implicit subgradient flow whose iterates can be determined explicitly.
The update with respect to u h , p h and s h can be regarded as proximal maps of E h,ε,δ (·, p j h , s j h ), E h,ε,δ (u j+1 h , ·, s j h ) and E h,ε,δ (u j+1 h , p j+1 h , ·), respectively. The following proposition guarantees stability and convergence of the algorithm.
holds. Particularly, Algorithm 2 terminates and E h,ε,δ (u j h , p j h , s j h ) j∈N is monotonically decreasing. Furthermore, every convergent subsequence Proof. The proof of the stability estimate follows as in the proof of [3, Thm. in the optimality conditions defining the iterates in Algorithm 2 and using the seperate convexity of E h,ε,δ . The stability estimate (4) implies that the sequence (u j h ) j is bounded. We choose a subsequence (u j h ) which converges to some u h,ε,δ ∈ S 1 (T h ). Since d t u j h → 0 it follows that u j +1 h → u h,ε,δ as well and by the equivalence of norms on finite-dimensional spaces we also have ∇u j h → ∇u h,ε,δ and ∇u j +1 h → ∇u h,ε,δ . Furthermore, since d t p j +1 h → 0 the optimality condition for the iterates (p j +1 h ) implies which means that p j h → (|∇u h,ε,δ | 2 + ε 2 ) 1/4 . The optimality condition for the iterate u j +1 h and d t u j +1 h → 0 then imply that for all v h ∈ S 1 (T h ). The optimality condition for the iterate s j +1 for all r h ∈ K ∩ S 1 (T h ). This means that for → ∞ we have for all r h ∈ K ∩ S 1 (T h ), i.e., (u h,ε,δ , s h,ε,δ ) satisfies the necessary and sufficient optimality conditions for a minimizer of E h,ε,δ .
Remark 8. Note that due to Proposition 4 all minimizers (u h,ε,δ , s h,ε,δ ) of E h,ε,δ are within the ball in L 2 (Ω) around the unique minimizer u h of E with radius c(ε + δh −2 ). Therefore, by Proposition 5 the iterates (u j h , s j h ) generated by Algorithm 2 get arbitrarily close to this ball around u h . 5.3. Regularization and splitting. Alternatively, we may consider only a regularization of the BV -seminorm and seek a minimizer u h,ε ∈ S 1 (T h ) of the functional In order to enforce the constraint strictly we can use an augmented Lagrangian ansatz, i.e., work with the regularized augmented Lagrangian functional to approximate the unique minimizer u h,ε ∈ S 1 (T h ) of E ε . We can use the following iteration to approximate the unique minimizer u h,ε .
The stability (and convergence) of Algorithm 3 follows from the theory of augmented Lagrangian methods (cf. [15,16]) using arguments from [22] which avoid differentiability assumptions on the involved functionals but exploit their (strong) convexity. (ii) In contrast to Algorithm 1 we do not include the term s j+1 h − s j h h in the stopping criterion since the optimality conditions we obtain in step (1) are coupled, i.e., as soon as η j+1 h − η j h = 0 the triple (u j+1 h , s j+1 h ; η j+1 h ) satisfies the optimality conditions for a saddle-point of L h,ε .
Since the minimizer in step (1) of Algorithm 3 is not directly accessible we consider for j ∈ N the augmented functional Noting that L h,ε is separately convex in u h , p h and s h we use as in the Heron-penalty method a semi-implicit subgradient flow to approximate a minimizer (u j+1 h , s j+1 h ) of L h,ε (u h , s h ; η j h ) in the (j + 1)-th iteration. Algorithm 4 (Decoupled regularized splitting method).
and then compute the minimizing (p j+1, +1 Otherwise, increase j → j + 1 and continue with (1). (2) of Algorithm 4 one can also stop the inner iteration using the stopping criterion

Remark 10. (i) Instead of the fixed number M of inner iterations in step
(ii) In our experiments we will set M = 1 and use the stopping criterion in step (4) of Algorithm 4. We will refer to this setting as the Heron-split method.
6. Application to piecewise constant segmentation. In this section we investigate the constrained minimization problem min u∈BV (Ω) with f ∈ L ∞ (Ω) and K := {v ∈ L 2 (Ω) : 0 ≤ v ≤ 1 a.e.}. This constrained minimization problem occurs, e.g., in the piecewise constant segmentation of images discussed in [13] where one aims at minimizing for a given image g ∈ L ∞ (Ω) the nonconvex functional with χ Σ being the characteristic function of the unknown set Σ. The authors propose an alternating minimization with respect to the unknown variables Σ, c 1 and c 2 . For a fixed set Σ the optimal c 1 and c 2 are given by However, the major difficulty amounts to solving the minimization problem defined by the functional M S(·, c 1 , c 2 ) for fixed c 1 and c 2 . Berkels et al. considered a relaxed functional of M S(·, c 1 , c 2 ) in [8], employed a finite element discretization and used the primal-dual method proposed by Chambolle and Pock in [12] using adaptive refinement techniques as well.
One can establish the existence of a minimizer u ∈ BV (Ω) ∩ K of E seg , however, the minimizer may not be unique. Restricting the minimization of E seg to S 1 (T h ) then also provides a minimizer u h ∈ S 1 (T h ). Analogously to Proposition 2 we obtain the following estimate. Lemma 6.1. Let u and u h be minimizers of E seg on BV (Ω) and on S 1 (T h ), respectively. We then have for d = 2 Proof. The proof is analogous to that of Proposition 2 and we therefore omit it.
The authors in [13] furthermore observed that the functional where θ(ξ) := max{0, 2|ξ − 1 2 | − 1}, has the same minimizers as E seg provided that β > α 2 f L ∞ (Ω) , cf. [13, Claim 1]. They employ a regularization of the BVseminorm and of the function θ and use an explicit gradient descent scheme to approximate a minimizer of E seg and E seg , respectively.
We remark that this procedure may be practical in this specific situation, however, the construction of the functional E seg only enables one to deal with this particular constrained minimization problem defined by E seg . Since we can adapt the iterative schemes presented in Section 5 to the minimization of E seg after making minor adjustments we will also present numerical experiments for the piecewise constant segmentation problem in the next section. We conclude this section with a result which is analogous to Proposition 4. Proposition 6. Let u h ∈ S 1 (T h ) and (u h,ε,δ , s h,ε,δ ) ∈ S 1 (T h ) × S 1 (T h ) be minimizers of E seg and of the functional respectively. Then we have Proof. Observing that for a ∈ R d we have |a| ≤ |a| ε ≤ |a| + ε we obtain which is the first inequality. Using Π K u h,ε,δ = s h,ε,δ and the arguments in the proof of Proposition 4 we obtain This proves the assertion.
Remark 12. Again, if we use a triangulation consisting of right triangles or tetrahedra with mutually orthogonal edges, we can relax the condition on δ as described in Remark 7 and may choose δ = h 1/2 /(α f ).
We first of all remark that for the chosen Ω, g and α the inclusion of the obstacle functions χ ≡ 1/5 and ψ ≡ 1/2 is not redundant. The solution for the classical ROF problem with these model parameters attains its maximal and minimal value which are greater than 3/4 and smaller than 1/10, respectively. Moreover, the scaling of the penalization parameter as δ = O(h) is justified by Remark 7 since the considered obstacles are globally constant and we use a triangulation consisting of halved squares. Particularly, Assumption 1 is satisfied. In Table 1 we displayed the iteration numbers needed by the algorithms to satisfy the stopping criteria (5), (6) and (7). Since the optimal choice of the step sizes σ 2 and σ for the split-split method and the Heron-split method, respectively, are not clear, we ran the algorithms for three different choices of step sizes. As one can observe the choice σ 2 = α and σ = α, respectively, lead to the smallest iteration numbers which is due to the balance of  fit-to-data and fit-to-constraint, i.e., when σ 2 = σ = α is chosen the split-split algorithm and the Heron-split method try to recover the datum g and to satisfy the bilateral constraint simultaneously. Note that the iteration numbers for σ 2 = 1 and σ = 1 for the split-split method and the Heron-split method, respectively, are astonishingly stable for different mesh sizes. However, this behavior cannot be inferred from the stability estimate given in Proposition 3. Regarding the Heronpenalty method, we ran the algorithm for the choices δ = h/α and δ = h. Obviously, due to a worse conditioning, the iteration numbers for the choice δ = h/α are much larger than for δ = h. However, for δ = h the iterates u j h generated by the Heronpenalty method stay far from the minimizer u h of E and the distance corresponds to the magnitude of α as predicted by Proposition 4. This effect is also illustrated in Figure 1 where the L 2 -error between the approximate minimizer u h of E generated by the split-split method and the iterates u j h and s j h , respectively, generated by the Heron-penalty method for both δ = h/α and δ = h is plotted against the number of iterations. For δ = h/α we see that the final iterates u N h and s N h serve as good approximations of u h and the approximation error is of order h. Yet, for δ = h, the iterates u j h do not have practical approximation properties and one would have to decrease the mesh size significantly to obtain a similar accuracy as for δ = h/α. Still, one can also observe that the different choices of δ do not considerably affect the approximation properties of the iterates s j h , i.e., one may set δ = h and choose the final iterate s N h as an approximation of u h . Yet, this is not justified theoretically by Proposition 4. Furthermore, the behavior of the active set The error curves for u h − u j h and u h − s j h for the split-split method are similar to those shown in Figure 1 for δ = h/α the only difference being that we have u h − u j h , u h − s j h → 0 as j → ∞ in contrast to the Heron-penalty method.
Example 2. The setting is as in Example 1 except for the parameter α, which is set to α = 500 here.
In Table 2 the iteration numbers for Example 2 are displayed. We see that the effect of higher iteration numbers for the split-split method and the Heron-split method for the choice σ 2 = σ = 1 and for the choice δ = h/α in the Heronpenalty method are now even more pronounced due to the choice of a large α. Moreover, Figure 2 shows once again that the iterates u j h generated by the Heronpenalty method for δ = h may not be used as a good approximation of u h since the constant α affects the L 2 -distance between u j h and u h while the iterates s j h seem to define accurate approximations of u h almost independently of the choice of δ. Again, the iteration numbers for σ 2 = 1 and σ = 1 appear to be stable.
The overall conclusion of the experiments is that the split-split method with σ 2 = α, Heron-penalty method with δ = h and the Heron-split method with σ = α lead to the smallest iteration numbers and have a similar performance. An advantage of the Heron-penalty method is that the choice of parameters is clear while the choice of the step sizes σ 2 and σ in the split-split method and the Heron-split method, respectively, remain unclear in general. However, the penalization parameter δ in the Heron-penalty method has to scale as h/α (or at least as h 1/2 /α to recover the convergence rate as stated in Remark 7) for the iterates u j h to be close to u h .  Table 2. Iteration numbers with (5), (6), (7) for Example 2 and σ 1 = h −3/2 for the split-split method and τ = 1 for the Heronpenalty and the Heron-split method.
7.2. Piecewise constant segmentation. In this section we report numerical experiments for the binary image segmentation model presented in Section 6. We consider the two-phase image segmentation of the image "cameraman". To this extent, we let g be the continuous piecewise affine function with coefficient vector whose entries consist of the (scaled) gray-values of the image ranging from 0 to 1. We further let α = 1500, χ ≡ 0, ψ ≡ 1, and consider the following algorithm of [13] which involves the solution of a constrained T V -minimization problem.
using either the split-split method or the Heron-penalty method or the Heron-split method.
Using a uniform triangulation consisting of halved squares with mesh size h ∼ 2 −8 we chose c 0 1 = 1 and c 0 2 = 0 and set the thresholding parameter to γ = 1/2. We set ε stop = 10 −3 in step (4). In step (2), we initialize the algorithms as in the preceding subsection and used ε = h 1/2 for the Heron-penalty and the Heron-split method. We set δ = h 1/2 in the Heron-penalty method and used σ 2 = σ = α f in the split-split and the Heron-split method and τ = 1 in the Heron-penalty and the Heron-split method. The stopping criteria in step (2) were chosen as follows: • Split-split method: • Heron-penalty method: • Heron-split method: Note that by multiplying the first inequality in (8) by h 3/2 the stopping criterion corresponds to ε (1) stop = h 3/2 10 −3 in Algorithm 1. We incorporated an h-dependent factor in ε (1) stop so that the gradient of u h is sufficiently well approximated. Figure 3 shows the original image and the outputs of Algorithm 5 using the splitsplit method, the Heron-penalty method and the Heron-split method in step (2), respectively, where we used the iterates s j h as approximations of a minimizer in step (2) in the Heron-penalty method. The three white horizontal lines in the outputs are due to image conversion. One should take into account that a uniform triangulation has been used based on the resolution of the original image (256×256) which makes the images seem to be slightly distorted. One can observe that all outputs are visually almost the same. The output of the split-split method and the Heron-split method visually do not differ at all, whereas the image generated by the Heron-penalty method has two black dots which do not appear in the other two images. These dots are not contained in the output of the Heron-penalty method when using δ = h 1/2 α f . The final values of c 1 and c 2 were c 1 = 0.5946 and c 2 = 0.1144 for both the split-split method and the Heron-split method and c 1 = 0.5951, c 2 = 0.1147 for the Heron-penalty method. The final values of the energy M S(Σ, c 1 , c 2 ) were 25.3745, 25.3863 and 25.3711 for the split-split method, the Heron-penalty method and the Heron-split method, respectively. The splitsplit method needed 5000 iterations in total for termination (6 outer iterations), while the Heron-penalty method needed 1382 iterations in total (7 outer iterations) and the Heron-split method needed 857 iterations in total (6 outer iterations). The Heron-penalty method with δ = h 1/2 α f needed 1675 iterations in total (6 outer iterations). Using u j h generated by the Heron-penalty method as approximations for the minimizer in step (2) led to the same results as for the case using s j h .   (2), respectively (horizontal white lines are due to image conversion).
8. Conclusion. We discussed a finite element discretization of a prototypical bilaterally constrained total variation minimization problem, proved convergence as h → 0 under a density assumption of C ∞ (Ω) ∩ K in BV (Ω) ∩ K, provided an experimental comparison between three approaches for the numerical solution of the discretized problem which differ in the numerical treatment of the nondifferentiable BV -seminorm and the nondifferentiable indicator functional encoding the constraint and addressed their stability. The numerical experiments reveal that the performance of the Heron-penalty method is negatively affected by the penalization parameter δ which determines how well u h,ε.δ approximates the discrete minimizer u h . However, we could have obtained lower iteration numbers if the penalization parameter had been chosen as δ = O(h 1/2 ) instead of δ = O(h). On the other hand, the split-split method avoids additional regularization of the nondifferentiable terms and yields low iteration numbers if the step sizes are chosen properly. The Heron-split method also yields low iteration numbers for appropriate step sizes, however, it lacks a stability estimate. Therefore, in our opinion, the split-split method seems to be the most appropriate out of these three algorithms for bilaterally constrained total variation minimization. Still, the choice of a proper step size in splitting methods is crucial and remains unclear. Here, it might be interesting to consider an adaptive step size adjustment to improve the convergence rate of splitting methods. The Heron-penalty method is also stable independent of the choice of the step size. The (optimal) scaling of the regularization and penalization parameters is clear but in most cases restrictive with regard to the penalization parameter δ which has to be chosen as δ = O(h 5/2 ) for d = 2 and δ = O(h 7/3 ) for d = 3. It will be interesting to investigate if this restrictive condition on the penalization parameter for the case of non-constant obstacles or triangulations with triangles which are not right-angled is really necessary. These topics as well as the application to total variation regularized damage evolution models is left for future research.