{'fraction_non_alphanumeric': 0.03963631938528936, 'fraction_numerical': 0.015419918988926168, 'mean_word_length': 4.049508426966292, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Phase field methods for binary recovery

We consider the inverse problem of recovering a binary function from blurred and noisy data. Such problems arise in many applications, for example image processing and optimal control of PDEs. Our formulation is based on the Mumford-Shah model, but with a phase field approximation to the perimeter regularisation. We use a double obstacle potential as well as a smooth double well potential. We introduce an iterative method for solving the problem, develop a suitable discretisation of this iterative method, and prove some convergence results. Numerical simulations are presented which illustrate the usefulness of the approach and the relative merits of the phase field models.


Introduction
A fundamental problem in the field of image processing is the following. Suppose we have a functionū defined on a bounded and piecewise smooth domain Ω ⊂ R N for N ≤ 3, which has been transformed by a linear operator S, and then corrupted by additive noise ζ, such that we have data The problem is to recoverū given y d . Two immediate issues are that (a) ζ is unknown, so we will not be able to findū even with a good model for the space in which it lies (b) inverting S may be ill-posed, so it will be difficult to find an approximation toū even if ζ = 0.
We investigate this problem in the case thatū is a binary function. We develop the theory with S an abstract operator, but in examples we take S to be the solution operator of an elliptic PDE. In this case the problem becomes one in PDE constrained optimal control.
Our approach to modelling the problem is to minimise an energy functional consisting of an L 2 fidelity term plus a phase field approximation to minimal perimeter regularisation. This can be thought of as a relaxation of the Mumford-Shah segmentation model. In our phase field approximation we use the Ginzburg-Landau functional with both the smooth double well and double obstacle potentials.

Motivating examples
First we give examples from both image processing and optimal control of PDEs motivating the study of this problem: • Image segmentation -We can represent a barcode by a 1D function which takes the value -1 when the barcode is white and 1 when it is black. When a barcode is scanned by a barcode reader this function becomes blurred (due to scattering in the air) and noisy (due to measurement error and imperfections in the barcode). So the machine only sees a corrupted signal, but from this it needs to determine the scanned barcode.
• Elliptic source recovery -Suppose we have noisy data of a quantity y, which is related to another quantityū by some physical law. For example, letū represent a heat source, then the long term temperature distribution y may be related toū by the solution of an elliptic PDE. Our goal could be to find the heat source that produces a particular temperature distribution.

Background material
For the above problems to be tractable we naturally require some knowledge of the form of the operator S and the noise ζ. We also usually assume a specific form ofū, as this influences the best model to use. For example, in the barcode problem we could assume that the function we are trying to recover is a binary function taking the values -1 and 1, and that the bars have a minimum width. Some sets of assumptions on S, ζ andū that are made in the literature are the following: 1. Denoising and deblurring -S is a blurring operator (maybe the identity), ζ is Gaussian noise, and u is a piecewise smooth function ( [ROF92], [CL97], [CE05]).
3. Binary image restoration -S is the identity, we have 'salt and pepper' noise, andū is a binary function ( [CEN06]). This kind of noise gives each point of a binary function a probability of switching to the other value, so the data y d is also binary.
Note that the above sets of assumptions have been named using terminology from image processing. Although our problem can be thought of as either an image processing or PDE constrained optimal control problem depending on the choice of S, we found most of the relevant literature to be from the image processing community. This is unsurprising since image processing is one of the main applications of binary recovery. We end up taking S to be the solution operator of an elliptic PDE, but try to use neutral language which reflects that our problem arises in these two fields.
For segmentation, which we focus on in this work, a large proportion of the literature modifies one of the following two models when formulating the problem of Section 1 mathematically. We now introduce these models so the reader can see how our approach fits with the existing literature.
• Model 1 (Mumford-Shah). This model, which was introduced in [MS89], looks for piecewise smooth functions that minimise an energy functional.
Let Ω i be disjoint open subsets with piecewise smooth boundaries such that the closure of Ω i is Ω. Let u be a function that is differentiable on Ω i , but which is allowed to be discontinuous across Γ := ∂Ω i \ ∂Ω. Then the Mumford-Shah model involves minimising where |Γ| denotes the N − 1 dimensional Hausdorff measure of Γ. The |Γ| term encourages minimising the length of the interface over which u is discontinuous.
If we restrict to minimising over binary functions that take the unknown value a i on Ω i (i = 0, 1), then this energy functional becomes For fixed Γ note that E 2 is minimised with respect to {a i } by setting So the problem reduces to just finding Γ, the locations of the discontinuities.
Due to the spaces of functions we are minimising over, both of the above variants of the Mumford-Shah model are nonconvex problems. In our work will use a relaxation of (1.1) based on a phase field approximation.
Sometimes the notation Ω |∇u| is used instead of |u| T V to highlight that the total variation of u is equal to this quantity when it is well defined. The first constraint in (1.2) says that the noise has zero mean and the second that it has standard deviation s.
BV (Ω, R) is the subspace of functions in L 1 (Ω) which have finite total variation. Minimising this model over u ∈ BV (Ω, R) can be related to the following problem for some value of σ: Note that (1.3) can be thought of as a relaxation of (1.1) with µ = 0; we minimise over a larger space of functions in order to get a convex problem.
If we restrict to minimising over binary functions then (1.3) becomes similar to the Mumford-Shah model. Suppose u only takes the known values a 0 < a 1 (i.e. u ∈ BV (Ω, {a 0 , a 1 })), then where the perimeter function Per(Σ) := Ω |∇χ Σ | and Γ is the set over which u is discontinuous. So for binary functions, total variation regularisation is equivalent to both perimeter regularisation and the interfacial length regularisation in the Mumford-Shah model. In fact (1.1) and (1.3) become equivalent.
Suppose that in addition to u ∈ BV (Ω, {a 0 , a 1 }) we have salt and pepper noise. Then the data is binary and both models reduce to the geometric problem min Σu⊂Ω |Σ u ∆Σ d | + σ(a 1 − a 0 )Per(Σ u ).
Here Σ u and Σ d denote respectively the sets where the unknown u and data y d take the value a 1 , |·| is now the N dimensional Hausdorff measure, and Σ u ∆Σ d is the symmetric difference between the sets.

Phase field model
We base our model on the Mumford-Shah model, but minimise over the space of functions BV (Ω, {a 0 , a 1 }), and generalise it to include the blurring operator S, which we suppose is known a priori. So we have the following nonconvex model with a parameter σ, which we will shortly relax in a different way to (1.3): We require S : L 2 (Ω) → L 2 (Ω) to be continuous, linear, and have the mean preservation property i.e. S1 = 1 and hence Sc = c for any constant function c. Continuity is helpful for proving existence of minimisers. Linearity and the mean preservation property allow us to recover a functionū : Ω → {a 0 , a 1 } from data y d by recovering a functionū : Ω → {−1, 1} from a scaled and shifted copy of y d , so long as a 0 and a 1 are known. We assume this to be the case and will therefore restrict our attention to a 0 = −1 and a 1 = 1 from now onwards.
Some examples of forms S could take are: 1. Solution operator of elliptic PDE -Let Su := y, where y solves the elliptic boundary value problem −α∆y + y = u in Ω ∂y ∂ν = 0 on ∂Ω. (1.5) For any u ∈ L 2 (Ω) this equation has a unique weak solution y ∈ H 1 (Ω) which satisfies the stability estimate where C s (α) := 1 1+α/Cp and C p is the Poincaré constant. So S has all the required properties. We also observe that evaluating S is well-posed, but inverting S is ill-posed, which motivates the need for our model. This is the operator we use for our numerics.

Convolution operator -Let
Su where φ α is a suitable probability distribution of 'size' α, for example the Gaussian distribution of mean zero and variance α, and * is the convolution operation. Such an operator is used in the barcode problem of [Ese04], [CG10] and [CVGO10].
In both of these examples we have a parameter α which controls the extent of the blurring effect. Large α corresponds to heavy blurring and small α corresponds to light blurring. In our work the value of α is known a priori since we assume complete knowledge of S. However there are applications where we may want to relax this assumption, for example the barcode problem of [Ese04]. In this application we do not know a-priori the distance of the barcode from the scanner, which means the level of blurring is unknown. This can be dealt with by fixing α to be some reasonable guess, or optimising for α at the same time as u.
We relax the model (1.4) by replacing the perimeter functional by the Ginzburg-Landau functional for some suitable Ψ : R → R, and then minimising over H 1 (Ω) instead of BV (Ω, {−1, 1}). So we consider arg min (1.7) We will focus on two different forms for the potential Ψ; the smooth double well potential and the double obstacle potential This approach, which is called a phase field approximation, results in a diffuse interface with minimisers no longer just taking the values {−1, 1}, but values in the interval [−1, 1]. It is still a nonconvex problem, but it has the advantage of allowing us to minimise over a smoother space of functions for which there is better developed theory. We are able to justify this approach with the following result.
A similar result holds for the double obstacle potential, and performing a calculation we get that c(Ψ 2 ) = π 2 (see [BE93]). To simplify notation we let σ i = σ/c(Ψ i ). This ensures that the weighting given to the regularisation is asymptotically σ for both potentials.
The different potentials lead to different formulations and we need to use different approaches to solve them. In particular, Ψ 1 leads to nonlinearity in the zeroth order terms, where as Ψ 2 causes nonlinearity by imposing constraints on the solution.

Literature review
We now mention other parts of the literature which overlap with aspects of this work.
Barcode problem. The 1D version of our problem is related to the barcode problem of Esedoglu in [Ese04]. This work was later extended by Choksi and Gennip in [CG10]. [CVGO10] uses similar ideas on QR barcodes. References for more general image processing literature can be found in Section 1.1.
PDE constrained inverse problems. A survey of the literature from the optimal control perspective can be found in [PS11]. In addition, [TC04] describes a number of applications where we want to recover piecewise constant functions, such as magnetic resonance imaging (MRI). The thesis [Hac06] discusses a wide range of techniques for geometric inverse problems. [TL07] recovers a piecewise constant diffusion coefficient from an elliptic PDE in 2D using the level set method.
Phase field. In [Hac06] there is a brief discussion of using a phase field approximation with the smooth double well potential for binary recovery. [Ese04] and [CG10] use this idea for numerical simulations, though they do not justify the approach analytically. Theory for the phase field approximation with the double obstacle potential can be found in papers by Blowey and Elliott, including [BE92] and [BE91b]. In [Sar10] the double obstacle potential is used in the context of image processing, but without deblurring.
Level set method. This is an alternative way of recovering the discontinuities in our problem. It is discussed in [TC04] and [TL07].
Approximation of Mumford-Shah. [CM99] and related papers prove Γ-convergence results for finite element approximations of the Mumford-Shah functional. These results have some relation to the convergence results that we obtain using a different approach.
Our work differs from existing work, and hence offers a new contribution, in the following respects: • We introduce the phase field approximation to the model right from the start (rather than at the last minute in order to allow numerical simulations). We therefore prove rigorous analytical results for this approximate model, which puts our approach on a much firmer footing than in existing work.
• Not only the smooth double well potential, but also the double obstacle potential is used for the phase field approximation. Results are proved for both simultaneously using an abstract framework.
• We thoroughly investigate the dependency of the model on the parameters and perform a systematic comparison of the smooth double well potential and the double obstacle potential on a 1D problem. This highlights some advantages and attractive features of the latter in this setting.

Layout
In Section 2 we introduce an abstract optimisation problem, an iterative method for finding critical points of this problem, and prove a convergence result for the iterative method. In Section 3 we show that (1.7) fits into this framework with both the smooth double well and double obstacle potentials. In Section 4 we discuss a gradient flow formulation of (1.7) and its link to the iterative method. In Section 5 we discretise the iterative method and prove another convergence result. We also look at a finite element discretisation for a particular choice of S. In Section 6 we demonstrate that implementations of the iterative method work well in 1 and 2 dimensions. In Section 7 the performance of using both potentials is compared in detail for a 1D problem. In Appendix A we describe how we choose the parameters in our model for the numerics.

Abstract framework
Rather than developing separate theory for solving (1.7) with the smooth double well and obstacle potentials, it is advantageous to introduce an abstract framework that both problems fit into. To this end let V and H be real Hilbert spaces with V compactly embedded in H, and let W be a closed convex nonempty subset of V . Let b : V × V → R and c : H × H → R be symmetric continuous bilinear forms with the properties Let l : V → R be a bounded linear functional and J : V → R a continuous convex functional. With these objects we can define the energy functional I : V → R by which for positive constants α 0 and C 0 we assume satisfies (2.1) Remark 2.1. The functional I can be decomposed in different ways into b, J, c and l.

Optimisation formulation
Consider the following optimisation problem: Find u ∈ W such that We can show existence of a solution to (2.2) with the following general result.
Proof. This follows from standard theory; we construct an infimising sequence which we know is bounded in V , so have a subsequence which weakly converges to an element of W , and this element is a minimiser of A by the properties of A 1 and A 2 .
Note that in general there is not a unique solution to (2.2).

Variational inequality formulation
By standard theory, solutions to (2.2) must satisfy the following: (2.3) Here we have used that J is a convex function, so it has a subdifferential ∂J, which by definition satisfies where ·, · denotes the duality pairing between V * and V . If J is in addition Gâteaux differentiable then (2.3) is equivalent to the following variational inequality: We often call solutions of (2.3) critical points of (2.2).
Remark 2.4. If c(η, η) ≤ κb(η, η) for all η ∈ V with κ < 1, then (2.3) has a unique solution. When we fit (1.7) into this framework, we find that this would require ε to be large. We intend to take ε small so that (1.7) approximates (1.4), which means we will not necessarily have uniqueness.
Note that solutions of (2.2) solve (2.3), but the converse is not necessarily true. We nevertheless aim to solve (2.3), as this is much easier in practice. Once a solution has been found, additional tests would have to be used to verify that the solution is a local minimiser of I.

Iterative method
We apply to (2.3) the following generalisation of the iterative method of Barrett and Elliott [BE91a]: (2.5) If J is in addition Gâteaux differentiable then this is equivalent to the following iterative method: Given (2.5) and (2.6) have unique solutions as they are equivalent to minimising a convex functional over W . Moreover we can prove the following convergence result.
Theorem 2.5. Every sequence {u n } generated by (2.5) satisfies and has a subsequence which converges in V to a critical point of (2.2) i.e. a solution of (2.3). Also, the limit of any subsequence of {u n } that converges weakly in V , and hence strongly in H, is a critical point of (2.2).
Proof. The proof is an extension to that of Theorem 6.1 in [BE91a]. To deduce (2.7) we test (2.5) with η = u n−1 and use the coercivity of b. Because of the assumptions on I, {u n } is uniformly bounded in V , so we can extract a subsequence, which we also denote by {u n }, that converges weakly in V and strongly in H to some element u ∈ W . The assumptions on b, c, l and J allow us to pass to the limit in (2.5) and deduce that u satisfies (2.3). The same argument applies to any subsequence, which proves the second part of the theorem. To see why the convergence in the first part of the theorem is strong in V , note that now we know u satisfies (2.3), we can combine this inequality with (2.5) to get The result then follows using the coercivity of b and the strong convergence of u n in H.

Binary recovery application
We now show that (1.7) with both the smooth double well and double obstacle potentials can be fitted into the framework of the previous section.

Smooth double well potential
Here and throughout this document (·, ·) denotes the L 2 (Ω) inner product. S * denotes the adjoint operator of S, which is defined as follows: For real Hilbert spaces U , V the adjoint operator of a continuous linear operator A : The above objects have the properties required in Section 2. Coercivity of b can be shown using a contradiction argument and that S0 = 0. J is well defined and continuous since H 1 (Ω) is continuously embedded in L 6 (Ω) for Ω ⊂ R N with N ≤ 3. I satisfies assumption (2.1) since and so Moreover I equals F ε from (1.7) with the smooth double well potential (up to an additive constant), so (2.2) becomes: Given y d ∈ L 2 (Ω) find arg min In this example we have an equality instead of a variational inequality because W is the full space V .
(2.6) gives the following iterative method for solving the above variational inequality, and it converges by Theorem 2.5: Given y d ∈ L 2 (Ω) and u 0 ∈ H 1 (Ω), for n = 1, 2, ... find u = u n ∈ H 1 (Ω) such that The above objects have the properties required in Section 2. As with the smooth double well potential, Moreover I equals F ε from (1.7) with the double obstacle potential (up to an additive constant), so (2.2) becomes: Given y d ∈ L 2 (Ω) find arg min Solutions to (3.3) satisfy (2.4), which becomes: Given y d ∈ L 2 (Ω), find u ∈ K such that (2.6) gives the following iterative method for solving the above variational inequality, which converges by Theorem 2.5: Given y d ∈ L 2 (Ω) and u 0 ∈ K, for n = 1, 2, ... find u = u n ∈ K such that (3.4)

Alternative iterative methods
In (3.2) and (3.4) the S * S term is taken implicitly, so we need to be able to invert the operator S * S −σ i ε∆ efficiently, otherwise these iterative methods will be too computationally expensive. In some cases this may be possible, for example if S is the identity, but in general this is not the case.
As we remarked earlier, the definitions of b and c that make I correspond to (3.1) and (3.3) are not unique. For example we can set b(u, η) = B(u, η) + ρ(u, η) and c(u, η) = C(u, η) + ρ(u, η) for some ρ ≥ 0. The ρ(u, η) terms cancel out in I, so defining B and C the same way b and c were defined earlier in this section gives the same optimisation problems (3.1) and (3.3). But the corresponding iterative methods are different. The point of this is that the ρ(u, η) term is convex (when η = u), so it gives us more flexibility in how we define B and C while still having b and c satisfy the coercivity and positivity assumptions.
In particular, for suitably large ρ we can take the S * S term explicitly (which in our framework corresponds to moving it from b to c), and also take the σi ε (u, η) term implicitly (i.e. move it from c to b). So for our examples this corresponds to taking A restriction such as ρ > max{ σi ε , C 2 s }, where C s is the stability constant from (1.6), is then sufficient for both b to be coercive and c to be nonnegative. So we have the following iterative methods, which are in general easier to solve computationally than (3.2) and (3.4).
Example 3.2 (Double obstacle). Given y d ∈ L 2 (Ω) and u 0 ∈ K, for n = 1, 2, ... find u = u n ∈ K such that When solving Example 3.1 in practice, it is more convenient for us to solve a linear equation. Therefore we linearise the J (u) term in (3.5) and consider the following iterative method.
This iterative method lies outside of our framework, so the convergence theory does not necessarily hold. However it works well in practice.
To finish this section we show how we can reformulate the iterative methods to remove S * (Su n−1 −y d ) when S is defined as in (1.5). For example, (3.6) becomes the following.
and y n−1 solves the weak form of (1.5) with u = u n−1 .

Gradient flow
In this section we investigate the gradient flow method for finding critical points of (3.1) and (3.3) from an initial guess u 0 . We prove that this method has some desirable properties, and note the link the to iterative method of the previous sections.

Smooth double well potential
Let u 0 denote the initial guess of the solution and consider the L 2 gradient flow of F 1 in (3.1).
Theorem 4.2. Problem 4.1 has a unique solution.
Proof. Note that Problem 4.1 is very similar to the Allen-Cahn equation with the smooth double well potential, and the proof follows using standard techniques. See for example the references in Theorem 4.5, where existence and uniqueness is proved for smooth potentials in order to show existence and uniqueness for the double obstacle potential in the limit.
Theorem 4.3. If u is a sufficiently smooth solution of Problem 4.1 then the energy F 1 (u(t)) decreases over time.
Proof. For some t ∈ (0, T ) we can test (4.1) with η = ∂ t u(t) to get Therefore as long as ∂ t u(t) is not zero almost everywhere we have and hence the energy decreases.

Double obstacle potential
We can formulate a gradient flow for F 2 from (3.3) in a similar way.
Problem 4.4. Given y d ∈ L 2 (Ω) and u 0 ∈ H 1 (Ω), find u ∈ K T with ∂ t u ∈ L 2 (0, T ; L 2 (Ω)) such that u(0) = u 0 and for all η ∈ K and almost all t ∈ (0, T ). Here Theorem 4.5. Problem 4.4 has a unique solution. Moreover, if u is a sufficiently smooth solution then the energy F 2 (u(t)) decreases over time.
Proof. This follows from a slight modification to the arguments for the double obstacle Allen-Cahn inequality in [BBG11, CE94, BGSS12, BE93, BE91b] to allow for the S * Su term.
For both potentials it is important to consider whether u(t) converges to a steady state as t → ∞. These types of issues are investigated in [Hal88], and in [CE94] for the 1D double obstacle potential. We do not discuss this as the focus of this work is on iterative methods. The scheme denoted by (2.5) of Barrett and Elliott motivated the convexity splitting implicit/explicit Euler scheme used in [ES93]. See also [Eyr98].

Discretisation
In this section we discretise the abstract iterative method of Section 2 in space and analyse convergence of the discretisation. We then apply this theory to a finite element discretisation of (3.5) and (3.6) for S defined by 1.5.

Discrete abstract framework
Suppose we have a family of subspaces V h ⊂ V and closed convex nonempty subsets W h ⊂ V h which approximate functions in W increasingly well as some parameter h → 0. In particular we suppose we have an approximation operator P h : W → W h such that and that every sequence Remark 5.1. Note that we do not require W h ⊂ W . If this holds then (5.2) follows automatically because W is a closed convex subset of a Banach space, and hence is weakly sequentially closed.
We now assume there exist objects b h , c h and l h which satisfy the same assumptions as b, c and l, with the boundedness and coercivity constants independent of h. We define and as in (2.1) we assume that there exist positive constants α 1 and C 1 independent of h such that So minimisers of I h over W h (which exist, since I h satisfies the same assumptions as I) satisfy the following discrete problem: If J is in addition Gâteaux differentiable then this is equivalent to the following discrete variational inequality: We need b h , c h and l h to approximate their continuous counterparts as h → 0. So we make the additional assumptions that for any bounded sequence {v h } ⊂ W we have as h → 0. With these assumptions solutions of the discrete variational inequality (5.4) approximate solutions of the continuous variational inequality (2.3) as h → 0, as the following theorem shows.
Theorem 5.2. For any sequence h n → 0 the sequence {u hn } of solutions to (5.4) has a subsequence which converges weakly in V , and hence strongly in H, to a critical point of (2.2) i.e. a solution of (2.3). Moreover, the limit of any subsequence of {u hn } that converges weakly in V , and hence strongly in H, is a critical point of (2.2).
Proof. For a given h we can find u h = arg min η h ∈W h I h (η h ), then for any η h ∈ W h , Here and throughout this section C denotes a generic constant independent of h which may vary from line to line. A similar result holds for l h , and c h is nonnegative, so By the boundedness of b and l, Combining this with (5.3) we get Now (5.1) and the continuity of J give that J(η h ) ≤ C. In addition (5.1) implies that for h less than some h 0 , η h V ≤ η V + C, and therefore u h V ≤ C.
From the above it follows that for any sequence h n → 0, {u hn } is bounded in V . So we can find a subsequence, which we also denote by {u hn }, that converges weakly in V and strongly in H to some u ∈ V . In fact u ∈ W by (5.2). We now show that u is a solution of (2.3).
Note that for all η ∈ W we have lim inf The final equality follows because lim n→∞ (b hn −b)(u hn , P hn η−u hn ) = 0 by (5.5) and lim n→∞ b(u hn , P hn η− η) = 0 by (5.1). The inequality follows from the lower semicontinuity of b(·, ·) and the continuity of b(·, η). Similar results hold for the c h and l h terms. This and the continuity and weak lower semicontinuity of J gives Hence u is indeed a solution of (2.3). The same argument applies to any weakly convergent subsequence, which proves the second part of the theorem.
Remark 5.3. We could also assume we have functionals J h satisfying the same assumptions as J, with the continuity independent of h, plus the additional property that v hn v in W for h n → 0 implies lim inf n→∞ J hn (v hn ) ≥ J(v). Then a proof almost identical to the above gives convergence for (5.4) with J replaced by J h . This allows numerical integration to be used on the J term.
As with (2.3) in Section 2, we can consider an iterative method for solving (5.4): Given u 0 h ∈ W h , for n = 1, 2, ... find u n h ∈ W h such that If J is in addition Gâteaux differentiable then this is equivalent to the following iterative method: Given u 0 h ∈ W h , for n = 1, 2, ... find u n h ∈ W h such that Since b h , c h and l h satisfy the same assumptions as b, c, and l, the above iterative method still has the energy decreasing property, and we get convergence of iterates to a solution of (5.4). Then as h → 0 the solutions of (5.4) converge to critical points of (2.2) by Theorem 5.2.

Finite element discretisation of (3.5) and (3.6)
Assume that Ω is polyhedral and let {T h } be a family of uniform regular triangulations of Ω into disjoint open simplices with a maximal element size h. Associated with each T h we have the piecewise linear finite element space where P 1 (T ) is the set of all linear affine functions on T . Also define K h := {v h ∈ V h : |v h | ≤ 1 in Ω} so that we have a finite element space analogous to K. Note that K h ⊂ K so Remark 5.1 applies. Take P h to be the operator that maps u ∈ W to the unique P h u ∈ W h such that This operator satisfies equation 5.1, see e.g. Chapter 2 in [Glo84].
Let S be the solution operator of (1.5), and denote by S h the discrete blurring operator. We intend this to approximate S, so we define S h to map u ∈ L 2 (Ω) to the unique y h ∈ V h satisfying (5.6) A stability estimate the same as (1.6) holds, so where as before C s (α) = 1 1+α/Cp . Also standard error analysis for elliptic PDEs says which combined with (5.7) gives that where S h is the discrete elliptic operator defined by (5.6), and y d,h is the L 2 -projection of y d onto V h . For ρ > max{ σ1 ε , C 2 s }, where C s is the stability constant from (5.7), all the assumptions of Theorem 2.5 are satisfied, so we get the decreasing energy property and convergence of iterates for the following discrete iterative method: Given y d,h , u 0 h ∈ V h , for n = 1, 2, ... find u h = u n h ∈ V h such that The assumptions of Theorem 5.2 are also satisfied, since for a weakly convergent sequence {v h } ∈ V we have Now using (5.7) and (5.8) we get and so (b h − b)(v h , ·) H 1 (Ω) * → 0 as h → 0 by the boundedness of v h V . Similar results hold for c h and l h . Therefore we have convergence of limit points of the above discrete iterative method to critical points of (3.1) as h → 0.
Remark 5.5. As mentioned before Example 3.3, when solving the smooth double well problem in practice, we solve a finite element discretisation of the linearised iterative method (3.7): Given y d,h , u 0 h ∈ V h , for n = 1, 2, ... find u h = u n h ∈ V h such that We use numerical integration on the linearised term. Note that the theorems do not necessarily hold for this iterative method, but it performs well in practice.
Example 5.6 (Double obstacle). Take the same definitions as in Example 3.2. In addition take V h as above, W h := K h , and define where S h is the discrete elliptic operator defined by (5.6), and y d,h is the L 2 -projection of y d onto V h .
For ρ > max{ σ2 ε , C 2 s } all the assumptions of Theorem 2.5 are satisfied, so we get the decreasing energy property and convergence of iterates for the following discrete iterative method: Given for all η h ∈ V h . Theorem 5.2 gives convergence of limit points of the above discrete iterative method to critical points of (3.3) as h → 0.

Algorithms
The discrete iterative methods in Examples 5.4 and 5.6 lead to the following algorithms for binary image recovery, which we implement and test in the next section.

Smooth double well potential
Given y d,h ∈ V h and an initial guess u 0 h ∈ V h , set n = 1 then: 1. Solve (5.10) for y n−1 h then p n−1 h ; 2. Solve (5.9) for u n h ; 3. If u n h − u n−1 h L 2 (Ω) < TOL terminate the algorithm. Else set n = n + 1 and go to step 1; An alternative stopping criterion would be to wait until the change in energy F 1 (u n h ) − F 1 (u n−1 h ) is sufficiently small. This has the advantage that the energy decreasing result then guarantees our algorithm terminates. However the stopping criterion in the above algorithm also gives a strong indication of a steady state, and it seems to work better in practice.
Note that despite the blurring and noise, y d,h still contains a lot of information about the solution. Therefore it makes sense to scale and threshold y d,h in order to get a good initial guess for u 0 h .

Double obstacle potential
The algorithm for this potential is the same as for the smooth double well potential, but we instead solve the variational inequality (5.11) in step 2. One method for solving the variational inequalities at each iteration is the primal-dual active set (PDAS) method. It is applied to solving the variational inequalities arising in the Allen-Cahn inequality in [BGSS12]. We implemented this method and found it to work well. However, for the numerics in the next sections we use an alternative method known as the Truncated Nonsmooth Newton Multigrid (TNNMG) method (see [Grä11,GK09]), which performs very well.

Numerics
In this section we show some numerical examples of binary recovery in 1 and 2 dimensions. The data is blurred by the solution operator of the elliptic PDE (1.5), with the parameter α controlling the level of blurring. It also has additive Gaussian noise of mean zero and variance γ.
We do the recovery using the discrete iterative methods of Remark 5.5 (based on the smooth double well potential) and (5.11) (based on the double obstacle potential). In practice we observe convergence of the full sequence of iterates to steady states, which are discrete critical points of (1.7). As we take ε and h small, we believe that these critical points closely approximate a global minimiser of the model (1.4). This is because the iterative methods give us discrete critical points of the approximate model (1.7), which seem to be at least discrete local minimisers of (1.7), as different initial iterates and (valid) values of ρ do not lead to different steady states. In addition, for small ε (and appropriate h) the critical points are close to being binary i.e. feasible minimisers of the model (1.4). We cannot be certain how close they really are to the global minimisers of (1.4) due to the lack of explicitly known global minimisers for interesting problems. Regardless, by artificially generating data from a known binary function, the numerical results show that for small ε (and appropriate h) our iterative methods are effective at recovering something close to the binary function.
The weighting given to the regularisation (the parameter σ), which defines the nonconvex model (1.4), is an important but challenging issue. If we take σ too small then recovered functions still have artifacts of the noise. If σ is too large then we loose some features we actually want to keep. We show some figures and discuss some results on the choice of σ for related problems in Appendix A.1, however the theory does not apply to our particular problem. In this section we just take values of σ that we have experimentally determined to work well for the problem at hand.
For the implementation we use the Distributed and Unified Numerics Environment (DUNE), see [BB07, BBD + 08b, BBD + 08a, BBD + 11, DKNO10,DKNO11]. DUNE provides interfaces for grids, solvers and finite element spaces. Therefore once the algorithms are implemented, it takes minimal effort to change features of the implementation that would usually be fixed, such as the grid type, the dimension of the problem, and the type of finite elements used.

1D numerics
The test problem in 1D is inspired by the barcode problem of [Ese04], which was mentioned as a motivating example in Section 1.1. We try to recover a binary function taking the values {−1, 1}, which one can imagine represents a cross section of a barcode (with values of -1 corresponding to black parts of the barcode and values of 1 corresponding to white parts). We suppose this binary function is corrupted, giving blurred and noisy data that we want to decode. The main difference between our test problem and the barcode problem in [Ese04] is that we have chosen blurring caused by the solution operator of an elliptic PDE instead of a convolution. Although this is not a realistic blurring operator specified by this application, if our approach is effective for this blurring operator then it is likely to be effective for other blurring operators.
The recovery using both the smooth double well and double obstacle potentials can be found in Figure 1. The black lines represent the binary function that we want to recover, the blue lines are the artificial data we generate by adding blurring and noise, and the red lines are the recovered functions for each potential. Even by eye it is not clear exactly how many 'bars' are in the binary functions, or the correct widths of the bars. But the recovered functions closely match the binary function we started with (which is why the black lines are almost hidden by the red lines), showing that our approach is effective. The figure also makes apparent one of the advantages of the double obstacle potential, which is that recovered functions take a form closer to what we actually want; binary functions.

2D numerics
The test problems in 2D involve recovering binary functions with discontinuities of various shapes. In this dimension the problems have a natural interpretation as deblurring and denoising of images, but we also view them as binary source recovery problems for elliptic PDEs. Figure 2 shows the recovery of a binary function using (1.7) with the smooth double well potential. The discontinuity is a 'blob' shape and is marked by a black line. The blurred and noisy data for this function is shown in Figures 2(a) and 2(b). Figure 2(c) shows the recovered function, with a yellow line marking the zero level set. We can see that the yellow line closely matches the black line, except for a slight mismatch at the concave parts of the discontinuity. Note that we cannot make the interface as small as for the 1D problem as the resolution of the grid needed to resolve it makes this computationally expensive. Our implementation is capable of adaptivity, which lessens this cost somewhat, but we will not demonstrate this functionality in this work. With this simple visualisation the recovered function using the double obstacle potential looks very similar, so we do not include a figure of it. Figure 3 (which can be interpreted in the same way as Figure 2) shows the recovery of a binary function with a letter 'A' shaped discontinuity. This time we use the double obstacle potential in (1.7), though the recovered function using the smooth double well potential looks similar. This example shows that the model can also recover discontinuities with corners reasonably accurately, but there is some rounding of these corners due to the regularisation.
To finish this section we show an example which relates to an application of binary image recovery in 2D. Figure 4 shows the recovery of a binary function representing a QR code with 25x25 blocks (the size typically used to encode a URL). The yellow lines mark the discontinuity of the binary function. Figure 4(a) shows the data with a red line marking the zero level set, and Figure 4(b) shows the recovered function. We see that features which are blurred below the zero level set (and which therefore would not be recovered by a simple projection) are nevertheless recovered by the model.

Comparison of potentials in 1D
Due to the Γ-convergence result of Theorem 1.1, we expect that critical points of (1.7) for a given value of σ using either the smooth double well or double obstacle potential will converge to critical point of (1.4) in the limit of small ε. Of course the critical points they converge to are not guaranteed to be the same, but agreement of the limits is observed in practice, and for very small ε the recovered functions for both potentials are almost indistinguishable. However it is well known that for phase field type problems, the interface should be well resolved in order for an accurate spatial approximation. This means that the smaller ε, the more grid points needed, and the higher the computational cost of the iterative methods. For many applications we only want to recover the location of the discontinuities in a binary function, which we suppose are given by the zero level set of the recovered function. This motivates us to consider in this section how well we can recover the locations of the discontinuities with ε of moderate size (rather than as small as possible), which is computationally cheaper. In this case the choice of potential does not just affect the implementation and speed of the iterative method; the recovered functions will in general look quite different, and there may be differences in how accurately or reliably the locations of the discontinuities are recovered. As in Section 6 we consider a problem with blurring caused by the solution operator of the elliptic PDE (1.5) and additive Gaussian noise of mean zero and variance γ. We use the discrete iterative method of Remark 5.5 for the smooth double well potential and (5.11) for the double obstacle potential.
At this stage it is helpful to recall the parameters we have introduced so far, as well as introduce a new parameter ω, the width of the smallest bar in the binary function. The parameters are contained in Table 7, and have been classified as follows: • Problem parameters -Define the problem we are trying to solve. In applications we have no control over these, though we suppose they are known a priori.
• Model parameters -Specify the model we will use to solve the problem. Different values can lead to the recovery of quite different functions, so they need to be chosen carefully.
• Approximation parameters -We do not work with the model, but rather an approximation of it. These parameters control how good the approximation is.
• Discretisation parameters -Affect the accuracy of the spatial discretisation in the iterative method.
• Iteration parameters -Determine the behaviour of the iterative method.
• Implementation parameters -Control the finer details of the implementation.
We also have a number of less significant implementation parameters that handle the imprecision of computer arithmetic. These will be set to sensible values and ignored in our discussion.
Motivated by the above discussion we now investigate differences between the smooth double well and double obstacle potentials in accuracy, reliability, speed, and implementational complexity.

Accuracy
Denote the binary function we want to recover byū and the recovered function by u ε,h . We measure the accuracy of the recovery by calculating the error quantity where P is the L 2 projection onto the space BV (Ω, {−1, 1}) (i.e. P (u ε,h ) = 1 when u ε,h ≥ 0 and −1 when u ε,h < 0). |u| T V is the total variation of u, as defined in Section 1.2. The integer part of E(u ε,h ) tells us the absolute difference between the number of bars in the projected recovered function andū. The decimal part tells us whether the discontinuities in the projected recovered function are in the correct locations. So E measures the accuracy of the recovery in a sense that matters in applications. We project because our best guess ofū should lie in BV (Ω, {−1, 1}). The downside of this is that P (u ε,h ) is not a minimiser of (1.7). It is important to note that the recovery using the double obstacle potential is naturally much closer to being binary than with the smooth double well potential, so projection is less necessary. This is a big advantage of using the double obstacle potential, which must be remembered when values of E(u ε,h ) seem comparable.
The test problems we use for our comparison use the same binary function as in Section 6.1 (which has ω = 1 113 ), and different levels of blurring and noise i.e. a range of values of α and γ. We first fix σ based on the size of ω (as described in Appendix A.1) then choose good values of the approximation and discretisation parameters (as described in Appendices A.2 and A.3). So we have σ = 1e − 4, ε = 7.06e − 4 and h = 2.77e − 4 for both potentials. Each realisation of the noise will be different, so we calculate an average E over multiple realisations of the noise. As we observed earlier, we get the same steady state of (1.7) regardless of the choice of iteration parameters. The same is true for implementation parameters. So we ignore both these types of parameters in our discussion of accuracy. We see in Figure 5 that neither potential is the most accurate in all circumstances. For moderate levels of noise (γ = 0.2), the double obstacle potential leads to a slightly more accurate recovery. However for high levels of noise (γ = 0.4), the smooth double well potential seems to perform slightly better. Without projection the double obstacle potential always leads to a recovery which is significantly more accurate than the smooth double well potential.

Reliability
By reliability we refer to the range of problems (i.e. the levels of blurring and noise) over which a binary function can be recovered with reasonable accuracy; as the amount of blurring and noise are increased, eventually the recovered function does not resemble the binary function we wanted. Note that this range will depend on σ. We do not do a detailed comparison of reliability, but feel that it is comparable for both potentials. For example, we can see in Figure 5 that α = 1e − 4 and γ = 0.4 is roughly the limit at which the correct number of bars can be recovered using either potential.

Speed
The time it takes to recover a function which resembles the binary function is an important practical consideration. Where as accuracy is independent of the implementation, this is certainly not the case for speed. All but the inner workings of each iterative method in our implementation are identical, so we will do our best to make a fair comparison of speed.
We perform this comparison for the binary function of Section 6.1, one choice of blurring and noise (α = 1e − 4 and γ = 0.2), and σ as in Sections 7.1 and 7.2. Choices of ε and h as well as iteration and implementation parameters have a big impact on speed, so we will test two different combinations of these parameters. Our timings can be found in Table 2.
The runtimes for 'accurate recovery' use ε and h as in Sections 7.1 and 7.2, and TOL as described in Appendix A.5. These values have been chosen to ensure robustness. The table also contains timings for 'rough recovery', where less conservative parameter values are used (ε = ω 2π , h = ω 20 , and TOL as described in Appendix A.5). For many problems we can still get a reasonable recovery with these parameter values, and it lowers the computation time significantly. The recovery times are comparable for each potential for both rough and accurate recovery, though the smooth double well potential has a slight advantage for this size of problem. However we remark that the recovery time of the double obstacle potential scales better as the number of degrees of freedom in the discretisation increases, so it has better performance in 2D.

Implementational complexity
Implementing the iterative method for the double obstacle potential is less standard as we are solving variational inequality rather than a PDE. But it is no more complicated than implementing adaptivity, which is needed for the computational cost of the iterative method for the smooth double well potential to scale well to dimensions 2 and higher.

Summary of comparison
Both potentials can accurately recover binary functions over the same range of blurring and noise. If no projection is used, the double obstacle potential produces significantly more accurate results. Even with projection it is more accurate for moderate levels of blurring and noise. Our implementation using the smooth double well potential is slightly quicker for both accurate and rough binary recovery on our 1D test problem. However our implementation using the double obstacle potential, which is overall no more complicated, scales better to many degrees of freedom and so tends to be quicker in higher dimensions.

A Parameter choices
In this appendix we describe our methodology for choosing parameter values for the numerical tests and comparisons in Sections 6 and 7.
A.1 Choice of model parameter σ We recover different functions for different values of σ, so it is important to choose the 'right' value. This is illustrated in Figure 6, where we show the recovered functions for the same problem as in Figure 1(a) for different values of σ. We see that σ = 5e − 3 leads to too few bars being recovered. The recovered function for σ = 1e − 6 follows the noise too much and does not resemble a binary function. With σ = 1e − 4 we recover something close to the binary function that generated the data, so we consider this to be a good value.
It is known that the choice of σ in (1.4) should be related to the variance of the noise. Noise with a large variance requires a large σ in order for good recovery. We could try to figure out the variance of the noise from the data and use this to choose σ, however there is not an explicit form for the relationship. Instead we choose σ based on the length scale of the features that we want to recover (i.e. the parameter ω), and use the same σ for all levels of noise. In applications this is generally known a priori e.g. for barcode recovery. This approach works well because we take σ be as large as possible while not removing the features we want to recover, and hence perform the maximum amount of denoising. We do not seem to pay a significant price for this large σ in cases where the noise is small, and this approach leads to a (a) σ = 5e − 3. simple rule for choosing σ. The literature that gives us a heuristic way of choosing such a σ is introduced below.
The following result shows that it is unwise to take σ too large.
But we also need to be careful not to take σ too small. In fact, since S is known we have the following result in the 1D case.
Another interesting result is Theorem 1.1 part 2 in [CG10], which proves more explicit conditions on σ to ensure exact recovery in the case that S is a convolution with a hat function in 1D. Due to our complicated form for S we are forced to use a more heuristic argument to choose a good value for σ. [CEN06] shows that for the 1D case in the absence of blurring and noise (i.e. binary data), local and global minimisers of (1.4) can be calculated explicitly for a given value of σ. These considerations suggest we should take σ to be smaller than a quarter of the size of the smallest object we want to recover. In particular, σ = ω 8 seems like a sensible choice. But this assumes binary data. We have blurring, which means the differences between the functions in the Su − y d 2 L 2 (Ω) term can be much smaller. Hence we take σ an order of magnitude smaller i.e. σ = ω 80 . This σ is still larger than the length scale of the noise (which is of order h), so the results in [CEN06] say it will be removed. Numerical experiments confirm that this choice of σ works well in practice.

A.2 Choice of ε
The phase field approximation in (1.7) results in solutions with interfaces of width o(ε). In order for an accurate spatial approximation we need a reasonable number of grid points across the interfaces. So a smaller ε requires more grid points and a higher computational cost. With this is mind we want to take ε as large as we can while still resolving the finest features of the binary function. So the choice of ε should be related to the value of ω.
We assume that there is a linear relationship between the optimal choice of ε and ω and deduce the constant of proportionality c 1 such that we get a good recovery with πε = c 1 ω. Note that πε is the asymptotic width of the interface for minimisers of the Ginzburg-Landau functional with the double obstacle potential, and a good approximation with the smooth double well potential. The width of interfaces in minimisers of (1.7), a perturbed Ginzburg-Landau functional, are approximately the same size. So c 1 can be thought of as the relative width of the interface compared to the width of the smallest bar.
To determine c 1 we recover a simple binary function which can be seen in Figure 7. We take ω 1 = ω 2 = ω 3 = 0.2 (i.e. bars of equal widths), as we found the case where all bars are at the finest length scale to be the hardest for accurate recovery. We consider different levels of blurring and noise and compute the error E of the recovered functions. We take σ to be the optimal value of ω 80 that we decided upon in Appendix A.1, and take πε = 50h to ensure that effects of the spatial discretisation do not distort our results.
We observe that for a high signal to noise ratio we can take c 1 very large and still get accurate recovery (α = 0.01 in Figure 8), even though the bars do not separate properly (see Figure 9(a)). For low signal to noise ratios (α = 0.1 in Figure 8) we need to take c 1 ≤ 0.5 for accurate recovery, though it is not until c ≤ 0.25 that the interfaces start to look reasonably sharp (see Figure 9(b)). As expected there is not an accuracy penalty for taking c 1 too small, however it increases computation time by forcing us to take smaller h in order to resolve the interfaces. This motivates us to take c 1 = 0.25 i.e. πε = ω 4 .
(b) c 1 = 0.25. Figure 9: The interfaces using the smooth double well potential with different values of c 1 . 9(b) shows the interfaces for c 1 = 0.25, which we decide is the optimal parameter value.

A.3 Choice of h
We use the same test problems as in Appendix A.2 to deduce a constant factor c 2 such that we get a good recovery with πε = c 2 h. Hence c 2 can be thought of as the number of grid elements across each interface. With a high signal to noise ratio (α = 0.01 in Figure 10) it can actually be advantageous to have few grid points across the interface. In this case the recovered function would have to deviate a long way from the binary function in order for the projection to take an incorrect value on even a single grid point, and the data does not force sufficient deviation. As a result we can actually get perfect recovery on coarse grids. However, if we have a poorly resolved interface we are not well approximating our model and we may get a bad recovery for low signal to noise ratios (α = 0.1 in Figure 10).
We do not want to adjust the relationship between ε and h for different levels of blurring and noise; we want a relationship for each potential that always works. This means we must properly resolve the interfaces. Figure 10 suggests that we can take c 2 = 5 for both potentials, however this leads to slightly jagged interfaces. Therefore we will again favour robustness and choose c 2 = 8 i.e. πε = 8h.

A.4 Choice of iterative parameter
The discrete iterative methods of Section 5.2 have valuesρ independent of h such that for all ρ >ρ the iterates decrease in energy and converge in some sense. For example, a possibleρ for the iterative method of Example 5.6 applied to the problem in Section 7.3 is max{ σ2 ε , C 2 s } = 0.999, where we use the Poincaré constant 1/π. However in practice we observe that the iterates of this method decrease in energy and converge for ρ ≥ 0.833. It is advantageous to take ρ small, as this results in fewer iterations and uses less total computational effort. So to maximise speed we experimentally determine a value of ρ which is as small as possible while still reliably giving a decrease in energy and convergence of iterates. This approach also works for the iterative method of Remark 5.5 for the double well potential, which lies outside of our framework. So for the speed comparison in Section 7.3 we use ρ = 0.833 for the smooth double well potential and ρ = 0.588 for the double obstacle potential. In the rest of the numerics, where speed is less of a concern, ρ is taken large (and larger thanρ if it is known) to ensure we get the expected behaviour of the iterative methods.

A.5 Choice of stopping criterion
We will never quite reach the steady state of the iterative method, so a decision needs to be made about when we are sufficiently close. For this purpose we use the stopping criterion introduced in Section 5.3 which terminates the algorithms when the L 2 norm of the difference between consecutive iterations is less than TOL.
Mostly we take TOL small so that we are effectively finding the exact steady state, but for the comparison of speed in Section 7.3 we need to avoid unnecessary iterations. Figure 11 suggests about 170 iterations will take us quite close to the steady state for the problem under consideration. This corresponds to taking TOL=3e−4 for the smooth double well and TOL=3.5e−4 for the double obstacle, and we use these values for the 'accurate recovery'.
In practice we just want a sufficiently accurate recovery as quickly as possible. Our feeling is that the binary function is usually sufficiently accurately recovered once the error is below 0.1. At this stage the correct number of bars have formed and the locations are probably known well enough (e.g. for a different algorithm to interpret the binary function as a barcode). We see in Figure 11 that the smooth double well potential achieves this in around 11 iteration, which corresponds to TOL=1.5e − 2. The double obstacle potential achieves this in around 9 iterations, which corresponds to TOL=4e − 2. We take these values for the 'rough recovery'.