AN ADAPTIVE FINITE ELEMENT METHOD IN L 2 -TV-BASED IMAGE DENOISING

. The ﬁrst order optimality system of a total variation regularization based variational model with L 2 -data-ﬁtting in image denoising ( L 2 -TV problem) can be expressed as an elliptic variational inequality of the second kind. For a ﬁnite element discretization of the variational inequality problem, an a posteriori error residual based error estimator is derived and its reliability and (partial) e  ciency are established. The results are applied to solve the L 2 -TV problem by means of the adaptive ﬁnite element method. The adaptive mesh reﬁnement relies on the newly derived a posteriori error estimator and on an additional local variance estimator to cope with noisy data. The numerical solution of the discrete problem on each level of reﬁnement is obtained by a superlinearly convergent algorithm based on Fenchel-duality and inexact semismooth Newton techniques and which is stable with respect to noise in the data. Numerical results justifying the advantage of adaptive ﬁnite elements solutions are presented.


Introduction
Image denoising is one of the fundamental tasks in mathematical image processing and aims at reconstructing an image from given noisy data. Often this is achieved by employing so-called variational methods which require to minimize an objective functional typically composed of a data fidelity and a regularization term; see, e.g., [44] for further details. Here we utilize an objective functional which was previously studied in [34] and several papers mentioned therein and is given by (P) minimize µ 2 krvk 2 L 2 (⌦) + 1 2 kv zk 2 L 2 (⌦) + ↵ R ⌦ |rv| l 2 dx over v 2 H 1 0 (⌦), where ⌦ ⇢ R 2 is the open, bounded and polygonal image domain, z 2 L 2 (⌦) denotes the given noisy image, and ↵ > 0 represents a regularization parameter. We always assume that 0 < µ ⌧ ↵ such that this model is a very close approximation to the renown total variation (TV) regularization model proposed by Rudin, Osher and Fatemi in their seminal work [35]: (1.1) minimize 1 2 kv zk 2 L 2 (⌦) + ↵ R ⌦ |Dv| over v 2 BV (⌦). Here, BV (⌦) denotes the space of functions of bounded variation and R ⌦ |Dv| the associated BVseminorm; see [22] for their definitions and further properties. For µ ! 0 in (P) it may indeed be shown that the associated sequence of solutions (u µ ) of (P) converges (weakly in L 2 (⌦)) to the solution u 0 of (1.1); see [34] and the references therein. Due to the importance of TV-regularization in edge preserving image restoration, many research e↵orts have been devoted to solving and extending the associated variational problem (1.1); see, e.g., [1, 12-14, 17, 33, 34, 37-40, 44] and the many references therein.
Due to given pixel frames, in order to solve the discrete version of (1.1) and/or its extensions or variations (like (P)) numerical solution schemes usually rely on finite di↵erence (FD) approximations where every image pixel corresponds to a degree of freedom in the discretization; see the monograph [44] and the references therein or, for instance, [12] or [34] for so-called primal-dual methods. In contrast to these uniform mesh approaches, in the present work we proposed to solve the variational problem resulting from the first order optimality characterization of the solution to (P) by means of an adaptive meshing technique.
Various discretization techniques including FD and finite element (FE) methods, respectively, admit adaptive techniques. Here we focus on adaptive finite element methods (AFEMs) due to their analytical advantages; see [3,6,43] for a general account of AFEM. Our interest in adaptive mesh refinement is motivated by the fact that it tries to achieve a given approximation quality of the numerical solution by ideally utilizing the smallest possible number of pixels. For instance, in our imaging context large homogeneous features, i.e. regions where u µ is (almost) constant, may be recovered well already by locally discretizing by a small number of triangles in a FE discretization of u µ . On the other hand, near edges and small scale features or texture the adaptive mechanism has to select a rather fine mesh automatically in order to guarantee a good approximation quality. In a posteriori residual-based AFEM this refinement process relies on local residual errors indicating triangles with large contribution to the overall error for refinement.
Returning to the above mentioned first order optimality conditions, setting where (·, ·) L 2 (⌦) denotes the usual L 2 (⌦)-inner product such that a : V ⇥ V ! R, with V := H 1 0 (⌦), becomes a V -coercive continuous symmetric bilinear form, l 2 V ⇤ , with V ⇤ = H 1 (⌦) denoting the (topological) dual space of V , and j : V !R := R [ {±1} is proper, convex and lower semicontinuous, the first order optimality characterization of the optimal solution u of (P) (where we keep µ fixed and thus neglect it in the notation for the solution) reads Note that in terms of the above notation the objective of (P) is given by where h·, ·i V ⇤ ,V denotes the duality pairing between V ⇤ and V . Since 1 2 kzk 2 L 2 (⌦) is a constant only, it may be neglected in the variational context such that (P) is equivalent to (1.4) minimize J(v) over v 2 H 1 0 (⌦). Observe that (1.3) represents a so-called variational inequality of the second kind [23]. Previously, in [9] a posterior residual based error estimates for a specific (and with respect to j a somewhat simpler) class of variational inequalities of the second kind were derived. In view of (1.3), our aim is to extend these estimates to situations where j models the TV-regularization term in image denoising and where the variational inequality is yet formulated in a di↵erent primal-dual fashion (see the subsequent section) which is in particular suitable for numerical solution schemes.
We mention that under a regularity assumption on the data, in [20] a priori error estimates for the FE solution of a regularized gradient flow for the TV-model (1.1) are studied. The adaptive finite element method is used in tomographic imaging such as di↵use optical absorption imaging [26] or fluorescence di↵use optical tomography [5,27] where either quadratic L 2 -or H 1 -regularization of Tikhonov-type (rather than the non-di↵erentiable TV-regularization) is applied. Some adaptive finite element approaches are used to solve a modification of the Perona-Malik model for image processing in [8] and [36]. Image segmentation based on adaptive finite elements is studied in [21], and e cient solvers relying on OcTree-based adaptive mesh refinement for parametric as well as non-parametric approaches in image registration can be found in [28][29][30].
The outline of the rest of the paper is as follows. Fenchel duality results useful for solving the optimization problem (1.4) are stated in Section 2. An a posteriori error bound to the finite element approximation is derived in Section 3. The e ciency of this bound is studied in Section 4. An adaptive finite element algorithm that utilizes this a posteriori error bound is proposed in Section 5. Numerical results of the adaptive finite element algorithm are presented in Section 6. Finally some conclusions are drawn.

Fenchel duality results
The optimization problem (P) can be solved by means of Fenchel duality; see [12,19,33,34]. In [34] it was shown that the Fenchel-dual problem of (P) is given by Here and below 'a.e.' stands for 'almost everywhere'. Letū be the solution of the problem (P) and letp be a solution of the dual problem (D). Then both quantities are related by the optimality system Note that the dual solutionp need not be unique in (2.1)-(2.2) which is due to the non-trivial kernel of the divergence operator occurring in weak form in (2.1). This fact may cause instabilities in primaldual or dual numerical solution schemes. In [34] this non-uniqueness is overcome by considering the regularized version of (P) given by and z 2 L 2 (⌦), and j(u) is replaced by the convex and Fréchet-di↵erentiable function j (u), which is defined as with > 0 a (small) regularization parameter. In [34] the dual problem of (P ) is shown to be The associated first order optimality condition relating the solutions of (P ) and (D ), respectively, is given by , for almost all x in ⌦; (2.5) see [34] and [24] for a related system in the context of Bingham fluids. We note that for small , 1 ↵ j is close to the BV-seminorm, Further, observe that j smoothes j locally around the kink of the norm-function involved in j. This local smoothing has a stabilizing e↵ect on primal-dual numerical solution schemes; see [34] for a related observation. Of course, other types of smoothing are possible without altering our subsequent considerations. An example, which is often found in the literature is q ✏ 2 + |rū| 2 l 2 . In the next section we derive residual-based a posteriori error estimators for a finite element approximation of the solution of the first order optimality system of either the original or regularized version of (P), respectively. By resorting to the formulation of the respective first order systems as a variational inequality of the second kind (in case of (P)) or a non-smooth variational system (in case of (P )), which can also be cast as a special variational inequality of second kind, both formulations will be treated within one framework. It, however, turns out that in case of (P ) stronger results are obtained, which is due to (2.5) (with the latter replacing the set-inclusion (2.2) by a equality relation).
For discretization purposes let V h ⇢ H 1 0 (⌦) =: V and P h ⇢ L 2 (⌦) =: Y be two finite dimensional subspaces defined by where ⇧ k denotes the set of polynomials with degree less than or equal to k defined on R 2 and T h denotes a regular triangulation of ⌦. For any open subset ! of ⌦ with Lipschitz boundary @!, we denote by L 2 (!), L 2 (!), H 1 (!) and L 2 (@!) the standard Lebesgue or Sobolev spaces equipped with the usual norms [2]. For any T 2 T h let E(T ) and N (T ) denote the set of its edges and vertices, respectively, and let us define Let h T be the diameter of T 2 T h , and let h E denote the length of E 2 E h . We denote by I h : V ! V h Clément's interpolation operator [16] satisfying the error estimates where the positive constants c 1 and c 2 depend only on the smallest angle in the triangulation; see [15,16,43].
The finite element approximationsū h 2 V h andp h 2 P h of the optimal solutionsū andp to (P) and (D) solve An analogous system is obtained when discretizing (2.4)-(2.5).

3.2.
A posteriori error estimate. Let e h :=ū ū h be the approximation error of the solutionū h of (3.5)-(3.6) to the true solutionū of problem (P). The bilinear form a : V ⇥ V ! R given by (1.2) is V -elliptic, i.e., there exist constants c a1 > 0 and c a2 > 0 such that Using the equivalence between the elliptic variational inequality of the second kind (1.3) and the optimization problem (P) and since V h ⇢ H 1 0 (⌦) we obtain Moreover, since F is convex the inequality holds true. Here we use fū h := f (ū h ). The convex function G need not be Fréchet di↵erentiable, but (in both cases considered above) it admits a generalized derivative which is given by In the case of G(·) = ↵ R ⌦ (·) dx becomes single-valued containing the Fréchet-derivative of j . Note that upon identification we have Y ⇤ = Y and h·, ·i Y ⇤ ,Y = (·, ·) L 2 (⌦) . Below, in a slight misuse of notation, we use gū h := g(rū h ). These derivatives together with (3.12) yield For our subsequent estimates we briefly recall the following relations.
, then it admits a well defined normal tracep ·ñ on = @⌦ which lies in H 1/2 ( ) [10,11,41] . Green's theorem implies whereñ @T denotes the unit outward normal to @T . Thus, from (3.16), Remark 1, (3.17) and the fact thatp h is constant on T , we infer Rearranging terms in (3.18), using the Cauchy-Schwarz inequality and applying the estimates (3.4), the following inequality is obtained: are defined on an edge E shared by the elements T and T 0 and represent the jump discontinuities in the approximation to the normal flux and the dual variable on the interface E, respectively; see, e.g., [3]. HereE(T ) andE h denote the interior edges sets, i.e.E(T ) : The inequality (3.19) together with (3.7) and (3.20) prove the next theorem.
Theorem 3. Letū 2 H 1 0 (⌦) be the solution of problem (P), letū h 2 V h be the finite element approximation to this solution withp h a corresponding solution of the dual problem (D), both given by the system (3.5)-(3.6). Then we have the following a posteriori error bound with some constant C > 0 depending only on c a1 and the shape regularity of the triangulation.
Based on our above findings we define the local error indicator . We end this section by a closer inspection of our estimates in the case where G(·) = R ⌦ (·) dx. In this case, the di↵erence kp h p k L 2 (⌦) is bounded above in terms of kū h ūk H 1 0 (⌦) as the next result shows.
) and its finite element approximation, respectively. Then, for any T 2 T h (⌦), we have Proof. We prove the result by partitioning ⌦ into subset according to where the max is taken in where we use |rū(x)|   |rū h | l 2 on M 2 and the second triangle inequality for the last estimate.
for some C > 0. Finally, recall that |p| l 2  ↵ and |p h | l 2  ↵ a.e. on T , respectively. Hence, where | · | denotes the two-dimensional Lebesgue measure. This concludes the proof. ⇤ From Theorem 3 and Lemma 4 we immediately conclude the following estimate.
Theorem 5. Letū 2 H 1 0 (⌦) be the solution of problem (P ), letū h 2 V h be the finite element approximation to this solution withp h a corresponding solution of the dual problem (D ). Then we have the following a posteriori error bound with some constant C > 0 depending only on c a1 and the shape regularity of the triangulation.

Efficiency of the a posteriori error estimator
The bound given in Theorem 3 is aimed to be used as a reference for an adaptive refinement process. In this section we show that the error bound is also reliable, i.e., up to data oscillations it also serves as a lower bound of the error.
We start by recalling a few results given in [43]. Let T i , i = 1, 2, 3, be the barycentric coordinates of T 2 T h (⌦) and b T := 27 Q 3 i=1 T i its associated triangle-bubble function. Let E i , i = 1, 2, be the barycentric coordinates of E 2 E h (⌦) and b E := 4 Q 2 i=1 E i its associated edge-bubble function. Then there exist constants c 3 , c 4 , c 5 , c 6 , c 7 that only depend on the shape regularity of the triangulation T h (⌦) such that for ⇢ T 2 ⇧ k (T ) and ⇢ E 2 ⇧ k (E), k 2 N 0 we have the following inequalities: where⇢ E is the extension of ⇢ E to the patchw E . For fixed E 0 2 E h (T )\{E} and for x 2 T we definẽ Let z h 2 V h be the finite element approximation to the given data z (see Remark 12) and let u h 2 V h andp h 2 P h be the finite element approximations to an optimal primal-dual pair (ū,p) which satisfy the system of equations (3.5)-(3.6). We proceed by stating and proving some preliminary results.
Lemma 6. There exists a positive constant C depending only on the shape regularity of the triangulation T h (⌦) such that Proof. For any v h 2 V h , according to the optimality conditions (2.4)-(2.5) and (3.5)-(3.6) we have and extend it by 0 to all of ⌦ (still keeping the notation f h for the extended function). Note that According to (4.2) and (4.3) we have the following inequalities: Using (4.10), (4.1) and the Cauchy-Schwarz inequality, we obtain the bound and consequently Finally, (4.8) follows from (4.13) and There exists a positive constant C depending only on the shape regularity of the triangulation T h (⌦) such that (4.14) Proof. We set ⇢ E := µ ⇥ @ūh @n ⇤ E + ⇥p h ·ñ ⇤ E and g h :=⇢ E b E , and we observe that b E 2 H 1 0 (! E ), which, like g h , can be extended by zero to a function in V h ⇢ H 1 0 (⌦) (again we keep the notation g h and b E for the extended functions). According to the optimality conditions (3.5)-(3.6) and Green's formula, we have Due to the regularity of the triangulation T h , the ratio h T /h E is uniformly bounded from below and above. Hence, (4.8) can be used in (4.16) to obtain (4.14). ⇤ Combining the previous auxiliary results, we conclude with the following theorem.
Theorem 8. Letū 2 H 1 0 (⌦) be the solution of the problem (P), letū h 2 V h the finite element approximation to this solution, and letp h be an approximation to a solution of the dual problem (D), given by the system of equations (3.5)-(3.6). Further let " T be given as in (3.23). Then there exists a positive constant C that depends only on the shape regularity of T h (⌦) such that Observe that in the case where G(·) = ↵ R ⌦ (·) Theorem 5 and Theorem 8 establish reliability and e ciency of the residual-based a posteriori error estimator for kū ū h k 2 for any > 0 (and, thus, for kū ū h k 2 H 1 0 (⌦) + kp p h k 2 L 2 (⌦) with a constant depending on 1 for estimating the di↵erence in the dual variables). In the case of G(p) = ↵ R ⌦ |p| l 2 dx, however, only an estimator for kū ū h k H 1 0 (⌦) is available where the e ciency bound depends on the di↵erence of the dual variablesp andp h . This is due to the non-uniqueness ofp in (2.1)-(2.2) and similarly forp h in the discrete analogue of this system. As a consequence, for decreasing h the quantityp h need not get closer top in the L 2 -norm. One only has that both vector-fields are bounded by ↵ in the pointwise almost everywhere sense and, hence, kp p h k L 2 (⌦)  2↵|⌦|.

Adaptive finite element algorithm
In the rest of the present work we consider the regularized problem (P ), whereū andp satisfy the optimality system (2.4)-(2.5). The finite dimension spaces V h ⇢ H 1 0 (⌦) and P h ⇢ L 2 (⌦), where the finite element solutionsū h andp h lie, are as follows These finite dimension spaces have the following bases The coe cients u i , i = 1, . . . , n, and p j , j = 1, . . . , 2m, are found by solving the system Note that (5.6) is non-smooth due to the presence of the max-term. For this reason we choose a generalized Newton method for solving (5.5)-(5.6). In fact, it turns out that the system (5.5)-(5.6) is semismooth and the generalized Newton iteration converges locally at a superlinear rate; see [34] for further details in the imaging context and [32] for a more general discussion of semismooth Newton methods.

5.1.
Adaptive finite element algorithm for image reconstruction. Throughout the remainder of this paper, the domain is given by ⌦ := (0, 1) 2 and the initial triangulation in our numerics is a uniform one consisting of 16 triangles.
Given a regular triangulation T h (k) and its corresponding finite element spaces V h (k) , P h (k) , we compute the finite element solutionsū h (k) 2 V h (k) andp h (k) 2 P h (k) by solving (5.5)-(5.6). We refer to [18,34] for the numerical solution algorithm which is of semismooth Newton type.
We observe that for such finite element spaces we obtain that ū h (k) = 0 and divp h (k) = 0, on each triangle. For each element T 2 T h (k) the local error indicator is given by Assumption 9. Based on the maximal resolution of the data image z, we assume that there exists a limit for the refinement, i.e. there exists h min > 0 such that for any T 2 T h (k) and all iterations k, we have h T > h min .

Bulk criterion.
For the local error for a triangle T 2 T h (k) given by " (k) T := ⌘ T + 1 2 ⌘ @T we observe that the orders of the element error ⌘ T and the jump discontinuity ⌘ @T are very di↵erent, respectively. The order of the element error is 10 7 , whereas the order of the jump discontinuity in the approximation to the normal flux is 10 and the one for the jump discontinuity in the approximation to the dual variable on the interface @T is 10 2 . Therefore we propose to use the bulk criterion for refinement; see [31]. At each iteration an element T 2 T h is marked to be divided if where, given the universal constants ✓ 1 > 0, ✓ 2 > 0 and ✓ 3 > 0, with ✓ i 2 (0, 1), i = 1, 2, 3, the sets M T h and M E h are such that (5.9) At each iteration, the resulting set of triangles T 2 T h (k) marked for refinement is denoted by D T h (k) . Remark 10. As it is well-known, the subdivision of an element T might lead to the division of one of its adjacent triangles as well in order to preserve the regularity of the triangulation. We refer to [42] for a description of possible mesh refinement techniques.

Noise level consideration.
The image data is given by z =û + ⌘ where ⌘ represents white Gaussian noise with mean zero and variance 2 . The solutionū for an appropriate regularization parameter ↵ satisfies the condition (5.10) kz ūk 2 L 2 (⌦)  2 , see [13]. Thus, if the variance 2 is known, an over-refinement in homogeneous image regions can be avoided by imposing the following constraint: T ,hmin) . Remark 11. The bound˜ 2 depends on the noise level 2 , on the size of the element T and on the resolution of the image given by h min . Numerically it needs to be computed in a similar way as it was done in Section 3.2 of [18]. For more details, we refer to Appendix A.

Final algorithm. Our proposed image restoration algorithm is as follows.
Algorithm AFEM-TV (1) Choose an initial regular mesh T h (0) .
(3) Compute, for each triangle T 2 T h (k) , the values ⌘ T and ⌘ E for all E 2E h (T ).

Remark 12.
The approximation of the data z given by z h is crucial for the results obtained by Algorithm AFEM-TV. In the case of noise free images, we use the following approximation where |T | is the area of the triangle T and (x c , y c ) are the coordinates of the center of gravity of T , as proposed in [4]. However in the case of noisy images, the information given by z(x c , y c ) is corrupted by noise and can no longer be used. Thus, instead of only considering z(x c , y c ) we consider z T = 1 |T | R zdx and we obtain Note that in our case we approximatez T by taking the average of the discrete data values of z contained in T .  kz T ,hmin) , which are subject to refinement. This increases the number of triangles marked for refinement in homogeneous regions. On the other hand, the result shown in the second row of Figure 2 relies on the condition (C ). It can be clearly observed that, although more triangles are used, the triangulation of the homogeneous region in the lower right corner remains coarser than in the first row. 5.5.2. Numerical specifications. In step (2) of algorithm AFEM-TV, the finite element solutions u h (k) 2 V h (k) andp h (k) 2 P h (k) are computed by means of the semismooth Newton algorithm described in [34]. The choice of the parameters for the semismooth Newton method is as follows: µ = 10 8 , = 10 6 , ↵ = 6 · 10 4 in the case of a noise free image and ↵ = 10 3 in the case of an image corrupted with white Gaussian noise with = 0.1.

6.
1. E ciency of the algorithm. We consider the noisy image related to the result in Figure 2.
The associated e ciency graph is presented in Figure 4. We observe that the algorithm using the adaptive mesh strategy outperforms the one using a uniform mesh. In order to reach an error of kū ū h k 2 L 2 (⌦) ⇡ 0.01, the adaptive mesh needs around 6000 triangles and the uniform mesh needs more than 130000; compare points PA4 and PU4. The numerical specifications are the ones given in section 5.5.2. The performance of this algorithm for more complex images is evaluated in the rest of the paper. 6.2. General remarks. The test images used in this section are the ones shown in Figure 5. The original images are the so called "Clock" in (a) and the "Cameraman" in (c), both of size 512 ⇥ 512 pixels. Thus, the maximal resolution of the image is h min = p 2/512. These images are corrupted by white Gaussian noise with standard deviation = 0.1; see (b) and (d). Unless otherwise indicated, the general numerical specifications in the following tests are as described in section 5.5.2. We present  Figures 7 and 9 we present the evolution of the peak signal-to-noise ratio (PSNR) according to the number of triangles.
6.3. Test of algorithm AFEM-TV. We start by reconstructing the original images of the Clock and the Cameraman presented in Figure 5 (a) and (c), respectively. The reconstruction of the Clock with Algorithm AFEM-TV is shown in the first row of Figure 6 (a)-(b) and the reconstruction of the Cameraman is in the second row (d)-(e). In both cases we observe that homogeneous regions like the background of the Clock or the sky in the Cameraman require only a coarse triangulation for their successful reconstruction. In fact, we observe that not only the sky of the Cameraman remains homogeneous but also that the shaded background of the Clock is well preserved. Detail regions like the portrait and the clock in the Clock or the camera in the Cameraman require more refinement. The edges also need more refinement, of course. In the resulting adapted meshes we can clearly observe the contours of the portrait, the clock and the book in the Clock and the contour of the man and the camera in the Cameraman. Small details like the numbers in the Clock are well recovered although the white buildings in the back of the Cameraman are not as well reconstructed. This is possibly due to the low contrast between the buildings and the sky. In Figure 6 (c), (f) we present the reconstructions given by the fine uniform mesh. We observe that although details like the buildings in the cameraman are better reconstructed with the uniform fine mesh, the cost of this reconstruction is extremely high when compared to the results of Algorithm AFEM-TV. In Figure 7 we show that a higher PSNR is obtained when using an adaptive mesh rather than using a uniform mesh with a similar number of degrees of freedom. Note here that AFEM is allowed to refine down to the scale given by the underlying pixel frame of the image. In the following experiment we restore the noisy images of the Clock and the Cameraman presented in Figure 5 (b) and (d). The results are depicted in Figure 8 and the performance of the algorithm is shown in Figure 9. The homogeneous regions like the background of the Clock or the sky in the Cameraman require a coarser mesh than the regions with smaller features, like the clock or the portrait in the Clock or the camera in the Cameraman. The edges are also well detected by the mesh which shows a stronger refinement at the contour of the book, the clock or the portrait in the Clock or the contour of the man in the Cameraman. Although the quality of the final reconstruction is not as good as the one given by the fine uniform mesh with h min,T = h max,T = 1/256, the reconstruction is still acceptable (compare, for instance, the PSNR-values shown in Figure 9). Similar to the reconstruction of noise free images, we observe in Figure 9 that a higher PSNR is obtained when using an adaptive mesh rather than when using a uniform mesh, both with approximately the same number of degrees of freedom.
6.4. Discussion on the convergence of the dual variablep h . Our numerics utilize the case where G = j in Section 3, i.e., the finite elements solutionsū h andp satisfy equation (5.6). Hence, according to Theorem 5, we observe that the a posteriori error bound also estimates the error in the dual variable given by kp p h k L 2 (⌦) since Lemma 4 shows that this error is bounded by the error of . Numerically it can indeed be shown that this error decreases as the number of triangles increases. As an example, in Figure 10 we show the convergence of the error kp p h k L 2 (⌦) for the reconstruction of the noise free and the noisy image of the Clock, respectively.

Conclusions
The results of this paper show that the adaptive finite element method (AFEM) may be used as an e cient discretization tool in total variation based image denoising. Based on a suitably chosen  Figure 9. Restoration of noisy images error estimator, it automatically adapts in image regions with details and edges while homogeneous regions are kept on a coarse discretization level. The proposed a posteriori error bound for adapting the mesh and finding a finite element approximation to the image is highly accurate for signals or images without noise. However, when noise is present, AFEM needs to be equipped with additional features to avoid over-refinement. For this purpose, condition (C ) of this paper relies on an estimate of the noise variance which takes the number of pixels (sample points) in each triangle into account and prevents an excessive refinement. A similar technique was used by the authors in earlier work for locally adapted regularization. Concerning the chosen finite element spaces we note that when µ = 0, then one has to resort to Raviart-Thomas FE spaces [10] as the pre-dual problem of (1.1) is Noise free clock then posed in H 0 (div) [33], the space of vector-fieldsp in L 2 (⌦) with divp 2 L 2 (⌦) andp ·ñ = 0 on = @⌦. Still, the refinement techniques and the tools for dealing with noise in the data etc. which are developed in the present paper will be important in this case. In this section we utilize the theory presented in [18] for finding local variance estimators.
We assume that the data z is a matrix of N ⇥ N pixels. Thus, in the discrete setting the noise ⌘ in z =û + ⌘ can be regarded as a sample of size N ⇥ N of a normally distributed random variable with zero mean and variance 2 . Each element T 2 T h (k) contains ! = d|T |/(h 2 min )e pixels. Then the random variable (A.1) W ! T = 1 2 P (i,j)2T (⌘ i,j ) 2 has the 2 -distribution with ! degrees of freedom, i.e. W ! T ⇠ 2 ! . Here, the notation (i, j) 2 T refers to the associated point of ⌦ lying in T . If u =û satisfies ⌘ = z û, then If, on the other hand, u is not the true signal and z u contains image details rather than noise only, then we expect that (A.3) S ! T = P (i,j)2T (z i,j u i,j ) 2 > P (i,j)2T (⌘ i,j ) 2 = 2 W ! T . Therefore, we search for a bound B such that S ! T > B implies that the element T contains details. Given N ⇥ N the total number of pixels in the image, we consider the expected maximum of a sample of size l = dN ⇥ N/!e coming from a 2 -distribution with ! degrees of freedom. We write E(max i=1,...,l W ! i ) =W +  . For a detailed description see [18] and [25].
Thus the bound˜ 2