Mixed Tikhonov regularization in Banach spaces based on domain ...
Applied Mathematics and Computation 218 (2012) 11583–11596
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Mixed Tikhonov regularization in Banach spaces based on domain decomposition Ivan Cimrák b, Valdemar Melicher a,⇑,1 a b
Department of Mathematical Analysis, NaM2 Research Group, Ghent University, Galglaan 2, 9000 Ghent, Belgium Department of Software Technologies, Faculty of Management Science and Informatics, University of Zˇilina, Univerzitná 8215/1, 01026 Zˇilina, Slovakia
a r t i c l e
i n f o
Keywords: Tikhonov regularization Domain decomposition Banach spaces Bounded variation Total variation
a b s t r a c t We analyze Tikhonov regularization where the forward operator is defined on a direct sum U of several Banach spaces U i ; i ¼ 1; . . . ; m. The regularization term on U is a sum of different regularizations on each U i . The theoretical framework, known for the case m ¼ 1, can be easily reformulated to the case m > 1. Under certain assumptions on weak topologies and on the forward operator it is possible to extend the results on the well-posedness of the minimization process, on its stability with respect to the data, and on the corresponding existence and convergence results, to the general case m P 1. We consider two particular regularizations, the total or the so-called bounded variation (BV) regularization and the smooth (H1 ) regularization. Assume that the domain X is decomposed into two non-overlapping subdomains X1 ; X2 . The direct sum of Banach spaces for m ¼ 2 will be U ¼ H1 ðX1 Þ BVðX2 Þ and the underlying regularization will be mixed H1 –BV regularization. We define proper weak topology on U and prove that the assumptions for the above mentioned framework are fulfilled. In such a way we get the theoretical results for a broad range of H1 –BV regularizations including the image denoising for the identity forward operator, the deconvolution problems for the convolution forward operator, or any other inverse problem for the continuous nonlinear forward operators. We demonstrate how the mixed regularization can be applied in practice. We suggest a denoising algorithm and compare its behavior with the BV regularization and the H1 regularization. We discuss the results of the denoising on several real world images justifying the usefulness of the mixed H1 –BV regularization. Ó 2012 Published by Elsevier Inc.
1. Introduction In many areas of inverse problems, the function to be reconstructed is expected to have different properties on different parts of the computational domain. As examples we mention magnetic resonance imaging or electrical impedance tomography where the image depicts different tissues with different physical and structural properties. Another example is image processing where the image contains smooth regions as well as sharp edges or textures. We use a domain decomposition technique that splits the domain into several parts and on each part the most suitable regularization is applied, chosen according to the expected level of smoothness of the solution. ⇑ Corresponding author. E-mail addresses: [email protected] (I. Cimrák), [email protected] (V. Melicher). Valdemar Melicher would like to acknowledge the support of GOA project BOF 07/GOA/006 of Ghent University, Ghent, Belgium and of BOF doctorassistant research mandate 01P09209T of Ghent University, Ghent, Belgium. 1
0096-3003/$ - see front matter Ó 2012 Published by Elsevier Inc. http://dx.doi.org/10.1016/j.amc.2012.05.042
11584
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
As an example of such an inverse problem let us consider image denoising problem where the forward operator is modeled by the identity. In the image restoration, two basic methods are widely used: The PDE based method and the energy method. We focus on the latter one. In this approach, an energy functional is constructed. This functional is minimized and the obtained minimum is the reconstructed image. Let us denote by X a domain in Rn ; n being the dimension. The following functional is to be minimized
T ðuÞ ¼ kf ukL2 ðXÞ þ RðuÞ:
ð1Þ
The functional typically includes the fidelity term kf ukL2 ðXÞ where f is the noisy data and u is an unknown function which is to be recovered. Further R is the regularization term, mostly it is a norm or a sum of several norms in some specific function spaces. The overview of the methods using different types of regularizations can be found e.g. in [3]. We briefly describe two choices of R, often forthcoming in the image processing problems. BV regularization. Rudin et al. in [18] suggested to use RðuÞ ¼ JðuÞ where JðuÞ is the seminorm in the space BVðXÞ of all functions with bounded variation, for the detailed definition see Section 4.2. Advantage: Since the functions from BVðXÞ have a small norm even if they have discontinuities along curves, this is a good regularization when the image is supposed to have sharp edges. Drawback: A serious drawback of this method is the staircase effect. The BV regularization favors constant functions with jumps over smooth functions. Therefore, wherever in the domain the solution should be smooth, the BV regularization introduces regions with constant intensities separated from each other by small and sharp edges. L regularization. For smooth functions, the L regularization can be used. In that case RðuÞ ¼ kLuk22 , where L is a differential operator, often Lu ¼ ru or Lu ¼ Du. In this case, RðuÞ is a seminorm in the usual Sobolev space Hk ðXÞ. The smoothing property of L can greatly remove the possible noise. Advantage: This type of regularization recovers smooth images very well. Drawback: The L regularization strongly penalizes discontinuous functions. Therefore the edges in images are smeared out. The result is oversmoothed and all the edges and the contours are lost. Mixed regularization. We propose the following approach. We decompose the domain into two non-overlapping subdomains such that one domain contains edges with their immediate neighborhood and the other subdomain contains the parts where the image is smooth. Then we apply the BV regularization on the first subdomain and the L regularization on the second. In such a way we avoid the drawback of the BV regularization in smooth parts and of the L regularization around the edges.The idea of using the BV regularization non-uniformly over the whole image has already appeared in the literature. We describe several methods, that use non-uniform regularizations. R Strong, Blomgren, Chan. For functions regular enough, the seminorm JðuÞ is equal to the following seminorm eJðuÞ ¼ X jruj. In [19] the authors first introduced a local adaptive algorithm for the BV regularization such that
T ðuÞ ¼ ~J x ðuÞ; with the constraint kf ukL2 ðXÞ ¼ r, where r is a noise level. The adaptive regularization term reads as
~J x ðuÞ ¼
Z
xðxÞjruj;
ð2Þ
X
where xðxÞ is a space variable weight function. This function is set inversely to the likelihood of the presence of an edge between any two neighboring discrete image locations. This allows for less regularization where edges are present and more regularization where there are no edges. This algorithm can be described as an adaptive BV regularization. Chan, Marquina, Mulet. Another approach presented in [6] leaves the BV regularization untouched and adds the space adaptive norm derived from an elliptic operator L resulting in the following functional to be minimized
T ðuÞ ¼ JðuÞ þ l
Z
UðjrujÞðLðuÞÞ2 ;
ð3Þ
X
where l is constant real number and U is a real function that indicates the presence of edges in the sense that its value approaches 0 when the gradient jruj is large, i.e., near edges. This algorithm can be described as a BV regularization with adaptive Hk regularization. Aubert, Vese. Next adaptive algorithm departs from the following functional
T ðuÞ ¼
Z
UðjrujÞ:
ð4Þ
X
There are many choices for U however they must conform to some hypotheses to give reasonable results in image recovery. These hypotheses have been formulated in [2] to satisfy the following principle in image analysis: The reconstructed image must be formed by homogeneous regions, separated by sharp edges. In [13] the authors defined similar hypotheses in order
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
11585
to conform to the more specific principle: The restored object can be decomposed into three scales – flat, gray and sharp edges. Charbonnier, in [7], presents many choices used in image reconstruction as well as a comparative study. Similar studies have been done in [5,13]. pffiffiffiffiffiffiffiffiffiffiffiffiffi The case when UðzÞ ¼ 1 þ z2 was studied in [10]. The minimization is considered under the constraint PR ¼ S where P R is the local variance and S is given a priori. This kind of algorithms can be described as an adaptive BV regularization. Gilles. The method presented in [11] is a generalization of the algorithm of Aujol, Aubert, Blanc-Féraud and Chambolle from [4]. The image is decomposed into three components u; v ; w representing structure, texture and noise, respectively. Noise is considered to be a rapidly oscillating function. The suitable function space for modeling of oscillating functions is the G space with the corresponding G norm, introduced by Meyer in [15]. It is known that the more oscillating function the lower G-norm this function has. Analysis of the minimization problems involving the combination of the BV and the G-norm regularization has been studied in [8]. So the texture belongs to the set Gl1 ðXÞ, where Gl1 ðXÞ :¼ u 2 GjkukG < l1 and the noise belongs to Gl2 ðXÞ where l1 l2 . The minimizing functional reads as
1 T ðuÞ ¼ JðuÞ þ kf u m1 v m2 wkL2 ðXÞ ; ð5Þ k with ðu; v ; wÞ 2 BVðXÞ Gl1 ðXÞ Gl2 ðXÞ. The functions m1 ; m2 are the weights with the following behavior: in a textured region, v must be favoured so m1 is close to 1 and m2 is close to 0. In an untextured region, it is opposite. This algorithm can be described as a BV regularization with adaptive two-stage G regularization.
In all above mentioned approaches, the regularization weights are positive over the whole X. In e.g. (5), the seminorm JðuÞ is evaluated over the whole domain. The decomposition of f to u and to the other components depends on the actual algorithm which is used to minimize the cost functional. Assume we have an a priori information on the image, that the BV component u should vanish in a part of the domain. Then it is preferable not merely rely on the algorithm to find the decomposition, but to a priori set the BV regularization constant to zero on this part.Therefore we decompose the domain of the image into several parts according to what regularity of the solution we expect in the particular part. Then we use different regularizations on the corresponding parts of the image. We locally decide which regularization is going to be used according to certain decision rules. Roughly speaking, if the reconstructed solution seems to have a jump, we use the BV regularization and otherwise we use the L regularization. 2. Contributions We develop a theoretical framework for the use of various regularizations in different parts of the domain. We employ a domain decomposition approach. This allows for zero weight of the specific regularization in the parts of the domain where this regularization is unwanted. The work is based on the results of [12]. These are extended to incorporate multiple regularizations. We assume that all the functional spaces involved are merely Banach spaces. Consequently it permits us to consider wider range of regularizations used in practice e.g. the BV regularization. Next, we employ the framework to image denoising by both the BV and the L regularization. Here the operator F modeling the forward problem is identity. However quite general nonlinear operators F are covered. The only restriction on F are those from Theorem 7. We show the well-posedness of the minimization problem and we verify the stability. Finally we show that a class of the minimization methods converges if the noise decreases to zero. Remark 1. One can object that our work strongly relies on the manner how the image is decomposed into the smooth part and the part containing the edges. But this dependence is common with all the algorithms involving the non-uniform regularizations. The function xðxÞ in (2), the function U in (3) and (4), the values l1 ; l2 in (5), these are all parameters that have to be chosen to tune the algorithms. Remark 2. Meyer [15] has introduced the space GðXÞ of functions that can be highly oscillating and still have small norm. This space is in some sense dual to the BV space. An image always have a component belonging to BVðXÞ and a component belonging to GðXÞ. The regularization RðuÞ ¼ kukGðXÞ is thus not suited for single use but always in the combination with the BV-norm. We aim at decomposing of an image into non-overlapping domains and on every domain we use single regularization. Therefore the G regularization does not fit our schema and we do not consider it.
3. Paper overview In [12], minimization of (1) is studied for a broad range of nonlinear operators, not only the identity. We can use these results to obtain the well-posedness and the stability of the minimization process. In Section 4 we extend the original assumptions from [12] regarding the operator and the function spaces setting to the case of multiple regularizations. At the end of the section we review the theoretical results.
11586
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
In Section 5 we verify that the above mentioned assumptions are valid for both regularizations we are interested in. We formulate then general result on the well-posedness and on the stability for a nonlinear operator with some specific properties. Further in Section 6 we discuss the details of the implementation and finally in Section 7 we present computational results on the benchmark images. 4. Theory The aim of this section is to formulate results on the well-posedness, on the stability and on the existence of the minimizing solutions to a general minimization problem analogical to (1). In Section 4.1 we rigorously define our problem for a general nonlinear forward operator. We introduce a functional analytical framework which covers the type of problems described in Introduction. In Theorem 1 we state theoretical results on the well-posedness of our problem, in Theorem 2 we verify the stability of the solutions on the noisy data. In Theorem 3 we state the existence of the minimizing solution and finally in Theorem 4 we get the convergence of the family of approximating solutions. In Section 4.2 we introduce the space BVðXÞ and we present some properties of this space. 4.1. General setting We take over the notations from [12]. Let us consider an operator equation modeling the forward problem
FðuÞ ¼ v ;
ð6Þ
where F : DðFÞ # U ! V is the corresponding forward operator, in general a nonlinear mapping between Banach spaces U and V. We study minimizers of the following functional
T a ðuÞ :¼ kFðuÞ v d krV þ aRðuÞ;
ð7Þ
þ
where r 2 R; 1 6 r < 1 and a 2 R is a global regularization parameter. We assume U can be decomposed into a direct sum of m Banach spaces U 1 ; . . . ; U m so that U ¼ U 1 U m , denoting by P i the corresponding projectors from U onto U i , i ¼ 1; . . . ; m. Assume that R is a linear combination of m regularizing functionals
RðuÞ ¼
m X bi Ri ðP i uÞ;
ð8Þ
i¼1
where Ri : U i ! ½0; þ1 is a convex and proper stabilizing functional with the domain
DðRi Þ :¼ fui 2 U i : Ri ðui Þ – þ 1g for i ¼ 1; . . . ; m. Real positive factors b1 ; . . . ; bm are assumed to be constant throughout this section. This setting is a generalization of that from [12], for m regularizing functionals. We make the following assumptions analogical to those from [12]. Assumption 1 (1) V and U i are Banach spaces such that there are topologies sV and sUi associated with these spaces, which are weaker than the norm topologies, for i ¼ 1; . . . ; m. Denote by sU the topology generated by sUi in the following way: uk ! u in sU topology when Pi uk ! Pi u in sUi topology, i ¼ 1; . . . ; m. Clearly, sU is weaker topology than the norm topology defined on the direct sum U. (2) k:kV is sequentially lower semi-continuous with respect to sV , i.e. if v k ! v with respect to the sV topology then
kv kV 6 lim inf kv k kV : k!1
(3) (4) (4) (6)
F : DðFÞ U ! V is continuous with respect to the topologies sU and sV . Ri : U i ! ½0; þ1 is proper, convex and sUi lower semi-continuous for i ¼ 1; . . . ; m. Note that also RðuÞ is then proper, convex and sU lower semi-continuous. DðFÞ is closed with respect to sU and D :¼ DðFÞ \ ðDðR1 Þ DðRm ÞÞ – ;. For every a 2 R; a > 0, arbitrary fixed non-negative values bi ; i ¼ 1; . . . ; m, and M > 0 the sets
Ma ðMÞ :¼ fu 2 D : T a ðuÞ 6 Mg
ð9Þ
are sU i sequentially compact in the following sense: every sequence fuk g in Ma ðMÞ has a subsequence, still denoted by fuk g such that fPi uk g converges in U i with respect to the sUi topology. Note that if this statement is true, we have then that every sequence fuk g in Ma ðMÞ has a subsequence, which is convergent in U with respect to the sU topology.
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
11587
From [12] we directly have all the theoretical results on the well-posedness and the stability of the solution to (7) and the convergence result of variational regularization methods implementing the minimization of (7) for the case when RðuÞ is of the form (8). Theorem 1 (Well-posedness). Assume that a > 0; there exists a minimizer of T a .
v d 2 V. Let F;
D; V and Ri ; U i for i ¼ 1; . . . ; m satisfy Assumption 1. Then
Theorem 2 (Stability). The minimizers of T a are stable with respect to the data v d . That means that if fv k g is a sequence converging to v d in V with respect to the norm topology, then every sequence fuk g satisfying
uk 2 argminfkFðuÞ v k krV þ aRðuÞ : u 2 Dg has a subsequence, which converges with respect to the ~ of T a as in (7). mizer u
ð10Þ
sU topology, and the limit of each sU convergent subsequence is a mini-
We recall the notion of a generalized solution in the Banach space setting. Definition 1. An element uy 2 D is called an R-minimizing solution if
Rðuy Þ ¼ minfRðuÞ : FðuÞ ¼ v g < 1: Note that the form of RðuÞ incorporates multiple regularizations. Therefore the properties of an R-minimizing solution strongly depend on b1 ; . . . ; bm . We discuss the proper choice of these parameters in Section 6. Theorem 3 (Existence). Let Assumption 1 be satisfied. If there exists a solution of (6) then there exists an R-minimizing solution. Let us denote by d the noise level in data v and by v d the corresponding noisy data, i.e. we assume kv v d kV 6 d. In the following theorem we characterize the conditions under which the minimizers of the cost functional (7) for noisy data v d converge to a minimizer of this functional when d ! 0. Theorem 4 (Convergence). Let F; D; V and Ri ; U i for i ¼ 1; . . . ; m satisfy Assumption 1. Assume that (1) there exists a solution to (6) (In [12] they assume the existence of a solution to (7) instead of solution to (6) but it is a mistake), (2) we have a sequence v k :¼ v dk 2 V strongly converging for dk ! 0 to v in V, (3) we have a decreasing sequence ak such that
ak ! 0 and
drk
ak
! 0;
Denote by uk a solution to (10) for a ¼ ak . Then a sequence fuk g has a convergent subsequence with respect to the sU topology and a limit of each sU convergent subsequence is an R-minimizing solution. Moreover if the R-minimizing solution uy is unique then uk ! uy with respect to sU . Further we recall some classical results from functional analysis. Lemma 1. The norms in Banach spaces are weakly lower semicontinuous. Theorem 5 (Kakutani). Let X be a Banach space. The following are equivalent: (i) X is reflexive (ii) The closed unit ball of X is weakly compact. (iii) Every bounded set admits a weakly convergent subsequence. Similar statement is due to [17]: Theorem 6. A normed linear space is reflexive if its unit ball is weakly compact. Further we rigorously define the BV space. Details about Hk ðXÞ can be found e.g. in [14]. 4.2. Space of BV functions 1 n n We will frequently make use of functions of bounded variation. We denote by C 1 c ðX; R Þ the set of functions in C ðR Þ with compact support in X.
11588
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Definition 2. By BVðXÞ we denote the subspace of functions u 2 L1 ðXÞ such that the following quantity is finite
Jðu; XÞ ¼ sup
Z X
n uðxÞdivnðxÞdx j n 2 C 1 c ðX; R Þ; knkL1 ðX;Rn Þ 6 1 :
We define the norm kukBVðXÞ ¼ kukL1 ðXÞ þ Jðu; XÞ. We point out that BVðXÞ endowed with the norm kukBVðXÞ is a Banach space and when X R2 we have that BVðXÞ L2 ðXÞ. From [1] we have the following results. Lemma 2. The norm k kBV ðXÞ is weakly lower semi-continuous with respect to the Lp topology for 1 6 p < 1. Lemma 3. Let S be a BV-bounded set of functions. Then S is relatively compact in Lp ðXÞ for 1 6 p < n=ðn 1Þ. S is bounded and thus relatively weakly compact in Lp ðXÞ for p ¼ n=ðn 1Þ and n P 2. 5. Domain decomposition method In the regularization theory, several regularizations are widely used to treat the ill-posedness of the underlying inverse problems. For the clear demonstration of our ideas, let us choose two particular regularizations: the L regularization and the total variation regularization, often called also the bounded variation (BV) regularization. Both regularizations have their drawbacks. The L regularization suffers from the oversmoothing when a discontinuous solution has to be recovered. The BV regularization appears to detect the discontinuities very well, however, when a smooth solution is to be recovered, a staircase effect appears. Therefore, if the unknown solution has discontinuities and on some parts of the domain is smooth, none of the two mentioned regularizations does recover this solution properly. As already mentioned in Introduction we follow a very simple idea: Let us use the L regularization on that part of the domain where we suppose the solution to be smooth and the BV regularization on that part of the domain where the discontinuities are likely to appear. In Section 4 we built up a theoretical framework for incorporation of several regularizations. The application we had on mind was the one of combination of several different regularizations on different parts of the computational domain. We rigorously define our inverse problem. Consider a domain X 2 Rn that is split into two non-overlapping subdomains X1 ; X2 so that X ¼ X1 [ X2 ; m ¼ 2. The first subdomain represents the part, where the solution is expected to be smooth, and the second subdomain should cover the discontinuities of the solution. At this point we do not exploit how this domain decomposition can be achieved. This topic is discussed in Section 6. The forward problem is represented by an (not necessarily linear) operator F operating between functional spaces U and V. For example, for the image recovery problems it is the identity operator, for the inverse problems in diffusion it can be an operator giving a solution of a partial differential equation. The only assumptions we put on F are those from Assumption 1. Further take U 1 ¼ Hk ðX1 Þ and U 2 ¼ BVðX2 Þ, where k is the desired order of expected regularity of the solution. We construct U by
U ¼ fu 2 Ln=ðn1Þ ðXÞ : ujX1 2 U 1 ; ujX2 2 U 2 g:
ð11Þ
Introducing the projector Pi : U ! U i being the restriction on Xi for i ¼ 1; 2 we have that U ¼ U 1 U 2 . We equip U with the following norm kukU ¼ kP1 ukHk ðX1 Þ þ kP2 ukBVðX2 Þ . Further take V ¼ L2 ðX1 Þ Lp ðX2 Þ, where p ¼ n=ðn 1Þ, the norm in V being kv kV ¼ kP 1 v kL2 ðX1 Þ þ kP2 v kLp ðX2 Þ . Consider the regularized least-square solution of the operator equation (6) with the noisy data v d for r ¼ 2
minfkFðuÞ v d k2V þ ab1 kP 1 ukHk ðX1 Þ þ ab2 kP2 ukBVðX2 Þ g: u2U
ð12Þ
Notice, that in the specific case F ¼ I, the minimization problem (12) is equivalent to two independent minimization problems: on X1 with Hk ðX1 Þ-regularization and on X2 with BVðX2 Þ-regularization. To analyse the minimization problem given by (12) we will apply the theoretical framework from Section 4. If the problem fulfills Assumption 1, all the theoretical results from Section 4 apply automatically. We verify Assumption 1 in the following subsection. 5.1. Verification of Assumption 1 (1) Consider U; U 1 ; U 2 as in (11). Let sV be topology generated by the weak norm on the direct sum L2 ðX1 Þ Lp ðX2 Þ. More specifically, the convergence of v k to v in sV means that P1 v k converges weakly to P 1 v in L2 ðX1 Þ and that P 2 v k converges weakly to P2 v in Lp ðX2 Þ. By sU1 denote the topology generated by the weak norm on Hk ðX1 Þ, by sU2 denote
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
the topology generated by the weak norm on Lp ðX2 Þ. So we have U and V the Banach spaces equipped with the topologies sU ; topologies. (2) Assume that v k ! v in V with respect to sV . By the definition it means that
P1 v k * P1 v P2 v k * P2 v
11589
sV which are weaker than the norm
in L2 ðX1 Þ; in Lp ðX2 Þ:
From Lemma 1, the strong norms in L2 ðX1 Þ and Lp ðX2 Þ are weakly lower semi-continuous. That results in
kP1 v kL2 ðX1 Þ 6 lim inf kP 1 v k kL2 ðX1 Þ ; kP2 v kLp ðX2 Þ 6 lim inf kP2 v k kLp ðX2 Þ : Putting it together we conclude
kv kV ¼ kP1 v kL2 ðX1 Þ þ kP2 v kLp ðX2 Þ 6 lim inf kP1 v k kL2 ðX1 Þ þ lim inf kP2 v k kLp ðX2 Þ 6 lim inf kv k kV ; which verifies that k kV is sequentially lower semi-continuous with respect to sV . (3) This assumption will be verified here only for the case of the identity operator F ¼ I with DðFÞ ¼ Hk ðX1 Þ BVðX2 Þ. We have to verify the continuity of the identity from DðFÞ to V ¼ L2 ðX1 Þ Lp ðX2 Þ with respect to weak topologies. Suppose uk ! u in sU . That means that
P1 uk * P 1 u in Hk ðX1 Þ; P2 uk * P 2 u in Lp ðX2 Þ: It directly follows that also
P1 uk * P1 u in L2 ðX1 Þ: This gives the continuity of the identity operator from DðFÞ to V with respect to the topologies sU and sV . If F is not the identity operator, we have to verify its continuity for each specific case separately. (4) Setting R1 ðuÞ ¼ kukHk ðX1 Þ and R2 ðuÞ ¼ kukBVðX2 Þ we directly have that RðuÞ ¼ b1 R1 ðP1 uÞ þ b2 R2 ðP2 uÞ is proper and convex. From Lemma 1 we know that R1 is sU1 lower semi-continuous. Lower semi-continuity of R2 with respect to sU 2 follows from Lemma 2. (5) Since DðFÞ is a direct sum of closed spaces, DðFÞ is closed too. (6) Suppose fuk g 2 Ma ðMÞ which means
kFðuk ÞkV þ kP 1 uk kHk ðX1 Þ þ kP2 uk kBVðX2 Þ 6 C: For reflexive space Hk ðX1 Þ we know from Theorem 5 that the unit ball is weakly compact. Since fP 1 uk g is bounded in Hk ðX1 Þ then there exists weakly convergent subsequence again denoted by fP 1 uk g. So we have that there exists u1 2 Hk ðX1 Þ such that P1 uk ! u1 with respect to sU1 topology. From Lemma 3 and the boundedness of fP 2 uk g in BVðX2 Þ we know that there exists subsequence again denoted by P2 uk convergent in sU2 topology with the limit u2 2 BVðX2 Þ. We construct u 2 U such that P 1 u ¼ u1 and P 2 u ¼ u2 and as a conclusion of the above preliminaries we state that from the sequence fuk g we are able to extract a subsequence convergent to u in sU topology. We successfully verified Assumption 1 for our setting. The parts (1), (2), (4)–(6) of Assumption 1 have been verified independently of F. The part (3) strongly depends on F and we verified it for the case of the identity operator F ¼ I. So we are able to formulate well-posedness and convergence results for general nonlinear F in the following theorem. Theorem 7. Consider the minimization problem (12). If F : DðFÞ U ! V is continuous with respect to the topologies sU and sV then Theorems 1–4 are valid for this minimization problem. Remark 3. It is possible to consider the BV semi-norm Jðu; X2 Þ instead of the BV norm kukBVðX2 Þ in all our considerations. The statements of Theorems 1–4 will remain valid after one more assumption on F, namely that F does not annihilate constant functions on X2 . The verification of the statements is analogous and will copy the reasoning in [1]. 6. Implementation In this section, we apply the ideas and the results obtained in Sections 4 and 5 to image denoising. The theoretical considerations from Section 5 are starting from the assumption that the image domain X is already split into two regions X1 ; X2 where X1 contains the smooth parts of the image and X2 the discontinuities. To determine the locations of the edges in the image, one can run any edge-detection algorithm that localizes steep variations of the image gray level. For an overview we
11590
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
refer to [20] and the references therein. Any of these techniques include an elaborated mechanism to choose the threshold according to which the decision is made whether edges are present or not at an image point. The lower the threshold, the more edges will be detected, and the result will be increasingly susceptible to noise, and also to picking out irrelevant features from the image. Conversely a high threshold may miss subtle edges, or result in fragmented edges. The aim of this paper is not to design an effective algorithm for edge detection but rather to design an algorithm that can denoise the image after we have the edges detected. Therefore we have focused our theoretical consideration on the denoising part having the image decomposed. For the actual decomposition we use the first order technique applied on the BV regularized image. We choose the BV-regularized image for the actual detection of the edges since the edges remain preserved after the BV regularization. Unlike in Sections 4 and 5 we use the BV seminorms Jðu; X2 Þ; Jðu; XÞ instead of the BV norms kukBVðX2 Þ ; kukBVðXÞ as justified in Remark 3. The forward operator is F ¼ I, the original image is v. The image v is polluted by noise of a known level d. The polluted image is denoted by v d . The ultimate goal of image denoising is to recover the original image v from v d . Obviously this is not possible exactly. We try to approximate v by a minimizer of the regularized cost functional (7). The quality of the reconstruction is determined by two factors: a and the ratio b1 =b2 . The parameter a in (7) determines the global weight of the regularization term RðuÞ. On the other hand, the parameters b1 ; b2 in (8) determine the relative weight of particular regularizations, in our case the Hk ðX1 Þ and the BVðX2 Þ regularizations. Consequently, the Hk regularization is weighted by c1 :¼ ab1 and the BV regularization is weighted by c2 :¼ ab2 . Given that c1 ; c2 are determined, Theorem 4 yields the convergence of a general method, satisfying the assumptions of the theorem, to the best-approximate solution in the sense of Definition 1. We minimize all the involved cost functionals by conjugate gradient method. We propose the following algorithm. The first and the second step comprise the decomposition of the domain and the actual denoising is done in the third step. Algorithm 1 1. Find u 2 BVðXÞ as a minimizer of
T BV ðuÞ :¼ kFðuÞ v d k2V þ cBV Jðu; XÞ:
ð13Þ
The regularization constant cBV will be determined by the discrepancy principle [16], see below. 2. Find the edges of u and decompose the image domain to X2 which contains all the edges and X1 :¼ X n X2 . 3. Find w 2 U as a minimizer of
T MIXED ðwÞ :¼ kFðwÞ v d k2V þ cMIXED kwkHk ðX1 Þ þ cBV Jðu; X2 Þ;
ð14Þ
where cMIXED is again determined via the discrepancy principle. 6.1. Discrepancy principle The discrepancy principle is a classical criterion used to determine regularization constant for ill-posed problems [9,16]. Solving FðuÞ ¼ v d is feasible only up to the level of noise, i.e we should not expect to find a regularized solution uda such that kFðuda Þ v d k < d (a is the corresponding regularization constant). In the same time one wants to introduce as much regularization as possible. When applying the discrepancy principle a is successively decreased untill kFðuda Þ v d k < Dd. The constant D 1 determines how much regularization is introduced and strongly influences the result. We however found the successive increase of a to be numerically more stable. Putting this together we propose the following simple update rule for a. 1. Set a to a small constant. 2. While kFðudak Þ v d k < Dd set akþ1 ¼ 2ak . We can get a more accurately by bisectioning on interval ak1 ; ak ½ for the final index k. 6.2. Domain decomposition To obtain the decomposition of the domain X into X1 and X2 from the step 2 of Algorithm 1 we employ the following simple approach. We compute the norm of the gradient of u. If it exceeds a certain value for a point of the image, this point belongs to X2 . Except this point we take into X2 also the points from some immediate neighborhood. The rest of image points belongs to X1 . 7. Numerical experiments We apply Algorithm 1 to eight real world 256256 images depicted in Fig. 1. Real images usually contain flat regions, smooth parts and textures. Usually none of these components vanishes at any point of an image. However, we assumed that we had an a priori information that some components are vanishing on a part of the domain. Namely we assume that the
11591
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Fig. 1. Images.
Table 1 Regularization constants and relative errors of L, of BV and of MIXED regularizations. Image
Caesar
Clown
Deer
Girl
Collection
Lighthouse
Lena
Airport
eL eBV eMIXED
0.0923 0.1128 0.0842 0.84 9.6 1.45
0.1047 0.1199 0.0876 0.675 8.4 2.657
0.1667 0.1411 0.1669 0.2225 1.63 0.888
0.0902 0.1094 0.0701 1.05 17.6 2.3133
0.1057 0.0904 0.1006 0.4275 3.22 0.7882
0.1009 0.1336 0.0952 1.6 29.28 2.87
0.1182 0.1272 0.0996 0.3975 3.8 3.7266
0.1085 0.1387 0.0997 1.39 24.96 6.6459
cL cBV cMIXED
images are composed of smooth regions separated from each other by sharp edges. This does not hold for general images. Consequently, one cannot expect that the domain decomposition of Algorithm 1 will be more powerful for general images than some algorithms from the Introduction. However, we will demonstrate that the domain decomposition works well. The discrepancy principle stopping constant D is set manually to DBV ¼ 0:5 for the BV regularization and DL ¼ DMIXED ¼ 1:3 for the smooth and the mixed regularizations (the regularization parameter cL of the smooth regularization is obtained in the same way as cBV by discrepancy principle). We use these values for all the experiments. First, we apply the uniform noise of 10% to all the images. In Table 1 one can find the corresponding numerical results: eL , eBV , and eMIXED are the relative L2 errors between the original image and the denoised one for the L, the BV and the MIXED regularization, respectively. cL ; cBV and cMIXED are the corresponding regularization parameters obtained by Algorithm 1. In majority of the cases the mixed regularization performs better than other regularizations and achieves the smallest relative error with respect to the original image. It is however misleading to compare the results on the basis of the relative errors only. Even if the reconstruction can be very close in L2 -sense, small details can be modified in a way that the image is visually degraded. For all the images, Algorithm 1 finds bigger cMIXED than cL even if DL ¼ DMIXED . It is caused by the following: If the L regularization is applied everywhere in the domain, it smears the edges. If cL was bigger than it is, the edges would be smeared much more and the reconstruction would be less accurate. Let us recall that discrepancy principle finds the maximal c’s such that the reconstructions are as close as possible to the noisy image, i.e. as one expects also as close as possible to the original one - up to the level of noise. When the L regularization is applied only on the smooth part of the domain X1 , it does not smear the edges in X2 . More L regularization can be applied in the smooth areas of the image and the reconstruction will be still close to the noisy image. The noise is then removed much better. Let us discuss three different images more closely. In Fig. 2 we present the L, the BV and the MIXED reconstructions for the Caesar image. The MIXED reconstruction is substantially better than the rest. First, the edges are recovered sharper than by the BV regularization (notice that the regularization constant cBV is the same). In the case of the MIXED regularization the BV
11592
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Fig. 2. Caesar: (a) original; (b) L-reconstruction; (c) BV-reconstruction; (d) MIXED-reconstruction.
Fig. 3. Caesar: (a) mask; (b) residuum of L-reconstruction; (c) residuum of BV-reconstruction; (d) residuum of MIXED-reconstruction.
regularization is applied only in some immediate neighborhood of the edges. Here the image is closely approximated by a function with a small BV norm - two flat regions divided by the edge. The gray scales are also recovered more accurately. The points in the smooth part of the domain X1 do not influence the intensities of points in X2 . The staircase effect is not visible. On the other hand, as already explained, more L regularization is applied in X1 . The noise is much better removed then in the case of plain L regularization. Fig. 3(a) displays the decomposition of X into X1 and X2 . White color represents the edges. i.e. X2 . Fig. 3(b)–(d) display the residua of the L, BV and MIXED reconstructions. The residua are normalized for clarity. The MIXED residuum contains clearly the least information, notice the details around Caesar’s nose and mouth.
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
11593
Fig. 4. Airport: (a) original; (b) L-reconstruction; (c) BV-reconstruction; (d) MIXED-reconstruction.
Fig. 5. Airport: (a) mask; (b) noisy image.
In Fig. 4 the results for the airport image are presented. This image contains a lot of small details. The MIXED reconstruction is clearly the best one. Even small details are recovered very well. Edges are sharp and intensities are correct. The corresponding decomposition of X and the noisy image are presented in Fig. 5.
11594
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Fig. 6. Lena: (a) original; (b) L-reconstruction; (c) BV-reconstruction; (d) MIXED-reconstruction.
Fig. 7. Lena: (a) mask; (b) residuum of L-reconstruction; (c) residuum of BV-reconstruction; (d) residuum of MIXED-reconstruction.
At last we present the results for Lena image. From Fig. 6 we can again conclude that the MIXED regularization performs best. The result seems a bit oversmoothed, e.g. the details of the feathers and on the face are lost. However the noise is removed much better then in the case of L regularization and the edges are sharper then for plain BV regularization. The corresponding domain decomposition and the residua are depicted in Fig. 7. The recovered images of Lena by MIXED regularization for different noise levels are depicted in Fig. 8. The respective numerical results are to find in Table 2.
11595
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Fig. 8. Lena: results obtained by mixed regularization for different noise levels of (a) 1%; (b) 3%; (c) 5%; (d) 20%.
Table 2 Lena: results obtained by mixed regularization for different noise levels of (a) 1%; (b) 3%; (c) 5%; (d) 20%. noise
1%
3%
5%
20%
eL eBV eMIXED
0.0158 0.0178 0.0138 0.0472 0.54 0.7786
0.0422 0.0481 0.0355 0.1838 2.06 9.1875
0.0635 0.0751 0.0541 0.3875 4.92 20.5375
0.1987 0.2642 0.1964 10.56 135.68 15.1287
cL cBV cMIXED
8. Conclusions We presented a theoretical framework for mixed Tikhonov regularization in Banach spaces. We first extended the results of [12], where the Tikhonov regularization in Banach spaces was tackled through completeness in a weaker than the norm topology. Particularly we looked for a solution in a direct sum of several Banach spaces and the weak topology was defined using projectors on these Banach spaces. Such projectors are naturally defined, if one considers a decomposition of the domain to certain subdomains and the Banach spaces are defined on these subdomains. The results obtained for domain decomposition based Tikhonov regularization are applicable for general possibly non-linear inverse problems, if the operator corresponding to the direct problem is continuous with respect to the weak topologies, see Theorem 7. We applied the framework to image denoising problem, when the direct problem is simply described by identity map. The image was decomposed to two parts. On one part, the smoothing L regularization was applied, while on the second part the BV regularization. We showed that the image denoising problem fulfills all the assumptions of the theoretical framework.
11596
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Finally, we implemented an algorithm which performs the mixed L–BV regularization. The results for a set of real world images confirm that the domain decomposition idea works well. Acknowledgment Large part of the presented work was done during the period when Ivan Cimrák was supported by the Fund for Scientific Research – Flanders FWO, Belgium. References [1] R. Acar, C.R. Vogel, Analysis of bounded variation penalty methods for ill-posed problems, Inverse Problems 10 (1994) 1217–1229. [2] G. Aubert, L. Vese, A variational method in image recovery, SIAM Journal of Numerical Analysis 34 (1997) 1948–1979. [3] J. Aujol, G. Gilboa, T. Chan, S. Osher, Structure-texture image decomposition–modeling, algorithms, and parameter selection, International Journal of Computer Vision 67 (2006) 111–136. [4] J.-F. Aujol, G. Aubert, L. Blanc-Féraud, A. Chambolle, Image decomposition into a bounded variation component and an oscillating component, Journal of Mathematical Imaging and Vision 22 (2005) 71–88. [5] P. Blomgren, T.F. Chan, P. Mulet, Extensions to total variation denoising, in: F.T. Luk (Ed.), Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, vol. 3162, 1997, pp. 367–375. [6] T. Chan, A. Marquina, P. Mulet, High-order total variation-based image restoration, SIAM Journal on Scientific Computing 22 (2000) 503–516. electronic. [7] P. Charbonnier, Reconstruction d’image: regularisation avec prise en compte des discontinuitis, Ph.D Thesis, University of Nice-Sophia Antipolis, 1994. [8] I. Cimrák, Analysis of the bounded variation and the G regularization for nonlinear inverse problems, Mathematical Methods in the Applied Sciences 33 (2010) 1102–1111. [9] H.W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Mathematics and its Applications, vol. 375, Kluwer Academic Publishers, Dordrecht, 1996. [10] G. Gilboa, N. Sochen, Y. Zeevi, Variational denoising of partly-textured images by spatially varying constraints, IEEE Transactions on Image Processing 15 (2006) 2281–2289. [11] J. Gilles, Noisy image decomposition: a new structure texture and noise model based on local adaptivity, Journal of Mathematical Imaging and Vision 28 (2007) 285–295. [12] B. Hofmann, B. Kaltenbacher, C. Pöschl, O. Scherzer, A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators, Inverse Problems 23 (2007) 987–1010. [13] K. Ito, K. Kunisch, Bv-type regularization methods for convoluted objects with edge flat and grey scales, Inverse Problems 16 (2000) 909–928. [14] A. Kufner, O. John, S. Fucˇík, Function spaces, Monographs and Textbooks on Mechanics of Solids and Fluids, Mechanics: Analysis, Noordhoff International Publishing, Leyden, 1977. Academia, Prague. [15] Y. Meyer, Oscillating patterns in image processing and in some nonlinear evolution equations, March 2001. The Fifteenth Dean Jacquelines B. Lewis Memorial Lectures. [16] V.A. Morozov, On the solution of functional equations by the method of regularization, Doklady Akademii Nauk SSSR 7 (1966) 414–417 (in Russian). [17] D.M. Oberlin, A measure-theoretic proof of a theorem on reflexivity, Proceedings of the American Mathematical Society 41 (1973) 325–326. [18] L. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1992) 259–268. [19] D.M. Strong, P. Blomgren, T.F. Chan, Spatially adaptive local feature-driven total variation minimizing image restoration, in: F. Preteux, J.L. Davidson, E.R. Dougherty (Eds.), Statistical and Stochastic Methods in Image Processing II, vol. 3167 of Proceedings of Society of Photo-Optical Instrumentation Engineers (SPIE), Bellingham, WA, 1997, Soc Photo Opt Instrumentation Engineers, SPIE – INT SOC OPTICAL ENGINEERING, pp. 222–233, Conference on Statistical and Stochastic Methods in Image Processing II, San Diego, CA, July 31– August 01, 1997. [20] D. Ziou, S. Tabbone, Edge detection techniques – an overview, International Journal of Pattern Recognition and Image Analysis 8 (1998) 537–559.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Mixed Tikhonov regularization in Banach spaces based on domain decomposition Ivan Cimrák b, Valdemar Melicher a,⇑,1 a b
Department of Mathematical Analysis, NaM2 Research Group, Ghent University, Galglaan 2, 9000 Ghent, Belgium Department of Software Technologies, Faculty of Management Science and Informatics, University of Zˇilina, Univerzitná 8215/1, 01026 Zˇilina, Slovakia
a r t i c l e
i n f o
Keywords: Tikhonov regularization Domain decomposition Banach spaces Bounded variation Total variation
a b s t r a c t We analyze Tikhonov regularization where the forward operator is defined on a direct sum U of several Banach spaces U i ; i ¼ 1; . . . ; m. The regularization term on U is a sum of different regularizations on each U i . The theoretical framework, known for the case m ¼ 1, can be easily reformulated to the case m > 1. Under certain assumptions on weak topologies and on the forward operator it is possible to extend the results on the well-posedness of the minimization process, on its stability with respect to the data, and on the corresponding existence and convergence results, to the general case m P 1. We consider two particular regularizations, the total or the so-called bounded variation (BV) regularization and the smooth (H1 ) regularization. Assume that the domain X is decomposed into two non-overlapping subdomains X1 ; X2 . The direct sum of Banach spaces for m ¼ 2 will be U ¼ H1 ðX1 Þ BVðX2 Þ and the underlying regularization will be mixed H1 –BV regularization. We define proper weak topology on U and prove that the assumptions for the above mentioned framework are fulfilled. In such a way we get the theoretical results for a broad range of H1 –BV regularizations including the image denoising for the identity forward operator, the deconvolution problems for the convolution forward operator, or any other inverse problem for the continuous nonlinear forward operators. We demonstrate how the mixed regularization can be applied in practice. We suggest a denoising algorithm and compare its behavior with the BV regularization and the H1 regularization. We discuss the results of the denoising on several real world images justifying the usefulness of the mixed H1 –BV regularization. Ó 2012 Published by Elsevier Inc.
1. Introduction In many areas of inverse problems, the function to be reconstructed is expected to have different properties on different parts of the computational domain. As examples we mention magnetic resonance imaging or electrical impedance tomography where the image depicts different tissues with different physical and structural properties. Another example is image processing where the image contains smooth regions as well as sharp edges or textures. We use a domain decomposition technique that splits the domain into several parts and on each part the most suitable regularization is applied, chosen according to the expected level of smoothness of the solution. ⇑ Corresponding author. E-mail addresses: [email protected] (I. Cimrák), [email protected] (V. Melicher). Valdemar Melicher would like to acknowledge the support of GOA project BOF 07/GOA/006 of Ghent University, Ghent, Belgium and of BOF doctorassistant research mandate 01P09209T of Ghent University, Ghent, Belgium. 1
0096-3003/$ - see front matter Ó 2012 Published by Elsevier Inc. http://dx.doi.org/10.1016/j.amc.2012.05.042
11584
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
As an example of such an inverse problem let us consider image denoising problem where the forward operator is modeled by the identity. In the image restoration, two basic methods are widely used: The PDE based method and the energy method. We focus on the latter one. In this approach, an energy functional is constructed. This functional is minimized and the obtained minimum is the reconstructed image. Let us denote by X a domain in Rn ; n being the dimension. The following functional is to be minimized
T ðuÞ ¼ kf ukL2 ðXÞ þ RðuÞ:
ð1Þ
The functional typically includes the fidelity term kf ukL2 ðXÞ where f is the noisy data and u is an unknown function which is to be recovered. Further R is the regularization term, mostly it is a norm or a sum of several norms in some specific function spaces. The overview of the methods using different types of regularizations can be found e.g. in [3]. We briefly describe two choices of R, often forthcoming in the image processing problems. BV regularization. Rudin et al. in [18] suggested to use RðuÞ ¼ JðuÞ where JðuÞ is the seminorm in the space BVðXÞ of all functions with bounded variation, for the detailed definition see Section 4.2. Advantage: Since the functions from BVðXÞ have a small norm even if they have discontinuities along curves, this is a good regularization when the image is supposed to have sharp edges. Drawback: A serious drawback of this method is the staircase effect. The BV regularization favors constant functions with jumps over smooth functions. Therefore, wherever in the domain the solution should be smooth, the BV regularization introduces regions with constant intensities separated from each other by small and sharp edges. L regularization. For smooth functions, the L regularization can be used. In that case RðuÞ ¼ kLuk22 , where L is a differential operator, often Lu ¼ ru or Lu ¼ Du. In this case, RðuÞ is a seminorm in the usual Sobolev space Hk ðXÞ. The smoothing property of L can greatly remove the possible noise. Advantage: This type of regularization recovers smooth images very well. Drawback: The L regularization strongly penalizes discontinuous functions. Therefore the edges in images are smeared out. The result is oversmoothed and all the edges and the contours are lost. Mixed regularization. We propose the following approach. We decompose the domain into two non-overlapping subdomains such that one domain contains edges with their immediate neighborhood and the other subdomain contains the parts where the image is smooth. Then we apply the BV regularization on the first subdomain and the L regularization on the second. In such a way we avoid the drawback of the BV regularization in smooth parts and of the L regularization around the edges.The idea of using the BV regularization non-uniformly over the whole image has already appeared in the literature. We describe several methods, that use non-uniform regularizations. R Strong, Blomgren, Chan. For functions regular enough, the seminorm JðuÞ is equal to the following seminorm eJðuÞ ¼ X jruj. In [19] the authors first introduced a local adaptive algorithm for the BV regularization such that
T ðuÞ ¼ ~J x ðuÞ; with the constraint kf ukL2 ðXÞ ¼ r, where r is a noise level. The adaptive regularization term reads as
~J x ðuÞ ¼
Z
xðxÞjruj;
ð2Þ
X
where xðxÞ is a space variable weight function. This function is set inversely to the likelihood of the presence of an edge between any two neighboring discrete image locations. This allows for less regularization where edges are present and more regularization where there are no edges. This algorithm can be described as an adaptive BV regularization. Chan, Marquina, Mulet. Another approach presented in [6] leaves the BV regularization untouched and adds the space adaptive norm derived from an elliptic operator L resulting in the following functional to be minimized
T ðuÞ ¼ JðuÞ þ l
Z
UðjrujÞðLðuÞÞ2 ;
ð3Þ
X
where l is constant real number and U is a real function that indicates the presence of edges in the sense that its value approaches 0 when the gradient jruj is large, i.e., near edges. This algorithm can be described as a BV regularization with adaptive Hk regularization. Aubert, Vese. Next adaptive algorithm departs from the following functional
T ðuÞ ¼
Z
UðjrujÞ:
ð4Þ
X
There are many choices for U however they must conform to some hypotheses to give reasonable results in image recovery. These hypotheses have been formulated in [2] to satisfy the following principle in image analysis: The reconstructed image must be formed by homogeneous regions, separated by sharp edges. In [13] the authors defined similar hypotheses in order
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
11585
to conform to the more specific principle: The restored object can be decomposed into three scales – flat, gray and sharp edges. Charbonnier, in [7], presents many choices used in image reconstruction as well as a comparative study. Similar studies have been done in [5,13]. pffiffiffiffiffiffiffiffiffiffiffiffiffi The case when UðzÞ ¼ 1 þ z2 was studied in [10]. The minimization is considered under the constraint PR ¼ S where P R is the local variance and S is given a priori. This kind of algorithms can be described as an adaptive BV regularization. Gilles. The method presented in [11] is a generalization of the algorithm of Aujol, Aubert, Blanc-Féraud and Chambolle from [4]. The image is decomposed into three components u; v ; w representing structure, texture and noise, respectively. Noise is considered to be a rapidly oscillating function. The suitable function space for modeling of oscillating functions is the G space with the corresponding G norm, introduced by Meyer in [15]. It is known that the more oscillating function the lower G-norm this function has. Analysis of the minimization problems involving the combination of the BV and the G-norm regularization has been studied in [8]. So the texture belongs to the set Gl1 ðXÞ, where Gl1 ðXÞ :¼ u 2 GjkukG < l1 and the noise belongs to Gl2 ðXÞ where l1 l2 . The minimizing functional reads as
1 T ðuÞ ¼ JðuÞ þ kf u m1 v m2 wkL2 ðXÞ ; ð5Þ k with ðu; v ; wÞ 2 BVðXÞ Gl1 ðXÞ Gl2 ðXÞ. The functions m1 ; m2 are the weights with the following behavior: in a textured region, v must be favoured so m1 is close to 1 and m2 is close to 0. In an untextured region, it is opposite. This algorithm can be described as a BV regularization with adaptive two-stage G regularization.
In all above mentioned approaches, the regularization weights are positive over the whole X. In e.g. (5), the seminorm JðuÞ is evaluated over the whole domain. The decomposition of f to u and to the other components depends on the actual algorithm which is used to minimize the cost functional. Assume we have an a priori information on the image, that the BV component u should vanish in a part of the domain. Then it is preferable not merely rely on the algorithm to find the decomposition, but to a priori set the BV regularization constant to zero on this part.Therefore we decompose the domain of the image into several parts according to what regularity of the solution we expect in the particular part. Then we use different regularizations on the corresponding parts of the image. We locally decide which regularization is going to be used according to certain decision rules. Roughly speaking, if the reconstructed solution seems to have a jump, we use the BV regularization and otherwise we use the L regularization. 2. Contributions We develop a theoretical framework for the use of various regularizations in different parts of the domain. We employ a domain decomposition approach. This allows for zero weight of the specific regularization in the parts of the domain where this regularization is unwanted. The work is based on the results of [12]. These are extended to incorporate multiple regularizations. We assume that all the functional spaces involved are merely Banach spaces. Consequently it permits us to consider wider range of regularizations used in practice e.g. the BV regularization. Next, we employ the framework to image denoising by both the BV and the L regularization. Here the operator F modeling the forward problem is identity. However quite general nonlinear operators F are covered. The only restriction on F are those from Theorem 7. We show the well-posedness of the minimization problem and we verify the stability. Finally we show that a class of the minimization methods converges if the noise decreases to zero. Remark 1. One can object that our work strongly relies on the manner how the image is decomposed into the smooth part and the part containing the edges. But this dependence is common with all the algorithms involving the non-uniform regularizations. The function xðxÞ in (2), the function U in (3) and (4), the values l1 ; l2 in (5), these are all parameters that have to be chosen to tune the algorithms. Remark 2. Meyer [15] has introduced the space GðXÞ of functions that can be highly oscillating and still have small norm. This space is in some sense dual to the BV space. An image always have a component belonging to BVðXÞ and a component belonging to GðXÞ. The regularization RðuÞ ¼ kukGðXÞ is thus not suited for single use but always in the combination with the BV-norm. We aim at decomposing of an image into non-overlapping domains and on every domain we use single regularization. Therefore the G regularization does not fit our schema and we do not consider it.
3. Paper overview In [12], minimization of (1) is studied for a broad range of nonlinear operators, not only the identity. We can use these results to obtain the well-posedness and the stability of the minimization process. In Section 4 we extend the original assumptions from [12] regarding the operator and the function spaces setting to the case of multiple regularizations. At the end of the section we review the theoretical results.
11586
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
In Section 5 we verify that the above mentioned assumptions are valid for both regularizations we are interested in. We formulate then general result on the well-posedness and on the stability for a nonlinear operator with some specific properties. Further in Section 6 we discuss the details of the implementation and finally in Section 7 we present computational results on the benchmark images. 4. Theory The aim of this section is to formulate results on the well-posedness, on the stability and on the existence of the minimizing solutions to a general minimization problem analogical to (1). In Section 4.1 we rigorously define our problem for a general nonlinear forward operator. We introduce a functional analytical framework which covers the type of problems described in Introduction. In Theorem 1 we state theoretical results on the well-posedness of our problem, in Theorem 2 we verify the stability of the solutions on the noisy data. In Theorem 3 we state the existence of the minimizing solution and finally in Theorem 4 we get the convergence of the family of approximating solutions. In Section 4.2 we introduce the space BVðXÞ and we present some properties of this space. 4.1. General setting We take over the notations from [12]. Let us consider an operator equation modeling the forward problem
FðuÞ ¼ v ;
ð6Þ
where F : DðFÞ # U ! V is the corresponding forward operator, in general a nonlinear mapping between Banach spaces U and V. We study minimizers of the following functional
T a ðuÞ :¼ kFðuÞ v d krV þ aRðuÞ;
ð7Þ
þ
where r 2 R; 1 6 r < 1 and a 2 R is a global regularization parameter. We assume U can be decomposed into a direct sum of m Banach spaces U 1 ; . . . ; U m so that U ¼ U 1 U m , denoting by P i the corresponding projectors from U onto U i , i ¼ 1; . . . ; m. Assume that R is a linear combination of m regularizing functionals
RðuÞ ¼
m X bi Ri ðP i uÞ;
ð8Þ
i¼1
where Ri : U i ! ½0; þ1 is a convex and proper stabilizing functional with the domain
DðRi Þ :¼ fui 2 U i : Ri ðui Þ – þ 1g for i ¼ 1; . . . ; m. Real positive factors b1 ; . . . ; bm are assumed to be constant throughout this section. This setting is a generalization of that from [12], for m regularizing functionals. We make the following assumptions analogical to those from [12]. Assumption 1 (1) V and U i are Banach spaces such that there are topologies sV and sUi associated with these spaces, which are weaker than the norm topologies, for i ¼ 1; . . . ; m. Denote by sU the topology generated by sUi in the following way: uk ! u in sU topology when Pi uk ! Pi u in sUi topology, i ¼ 1; . . . ; m. Clearly, sU is weaker topology than the norm topology defined on the direct sum U. (2) k:kV is sequentially lower semi-continuous with respect to sV , i.e. if v k ! v with respect to the sV topology then
kv kV 6 lim inf kv k kV : k!1
(3) (4) (4) (6)
F : DðFÞ U ! V is continuous with respect to the topologies sU and sV . Ri : U i ! ½0; þ1 is proper, convex and sUi lower semi-continuous for i ¼ 1; . . . ; m. Note that also RðuÞ is then proper, convex and sU lower semi-continuous. DðFÞ is closed with respect to sU and D :¼ DðFÞ \ ðDðR1 Þ DðRm ÞÞ – ;. For every a 2 R; a > 0, arbitrary fixed non-negative values bi ; i ¼ 1; . . . ; m, and M > 0 the sets
Ma ðMÞ :¼ fu 2 D : T a ðuÞ 6 Mg
ð9Þ
are sU i sequentially compact in the following sense: every sequence fuk g in Ma ðMÞ has a subsequence, still denoted by fuk g such that fPi uk g converges in U i with respect to the sUi topology. Note that if this statement is true, we have then that every sequence fuk g in Ma ðMÞ has a subsequence, which is convergent in U with respect to the sU topology.
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
11587
From [12] we directly have all the theoretical results on the well-posedness and the stability of the solution to (7) and the convergence result of variational regularization methods implementing the minimization of (7) for the case when RðuÞ is of the form (8). Theorem 1 (Well-posedness). Assume that a > 0; there exists a minimizer of T a .
v d 2 V. Let F;
D; V and Ri ; U i for i ¼ 1; . . . ; m satisfy Assumption 1. Then
Theorem 2 (Stability). The minimizers of T a are stable with respect to the data v d . That means that if fv k g is a sequence converging to v d in V with respect to the norm topology, then every sequence fuk g satisfying
uk 2 argminfkFðuÞ v k krV þ aRðuÞ : u 2 Dg has a subsequence, which converges with respect to the ~ of T a as in (7). mizer u
ð10Þ
sU topology, and the limit of each sU convergent subsequence is a mini-
We recall the notion of a generalized solution in the Banach space setting. Definition 1. An element uy 2 D is called an R-minimizing solution if
Rðuy Þ ¼ minfRðuÞ : FðuÞ ¼ v g < 1: Note that the form of RðuÞ incorporates multiple regularizations. Therefore the properties of an R-minimizing solution strongly depend on b1 ; . . . ; bm . We discuss the proper choice of these parameters in Section 6. Theorem 3 (Existence). Let Assumption 1 be satisfied. If there exists a solution of (6) then there exists an R-minimizing solution. Let us denote by d the noise level in data v and by v d the corresponding noisy data, i.e. we assume kv v d kV 6 d. In the following theorem we characterize the conditions under which the minimizers of the cost functional (7) for noisy data v d converge to a minimizer of this functional when d ! 0. Theorem 4 (Convergence). Let F; D; V and Ri ; U i for i ¼ 1; . . . ; m satisfy Assumption 1. Assume that (1) there exists a solution to (6) (In [12] they assume the existence of a solution to (7) instead of solution to (6) but it is a mistake), (2) we have a sequence v k :¼ v dk 2 V strongly converging for dk ! 0 to v in V, (3) we have a decreasing sequence ak such that
ak ! 0 and
drk
ak
! 0;
Denote by uk a solution to (10) for a ¼ ak . Then a sequence fuk g has a convergent subsequence with respect to the sU topology and a limit of each sU convergent subsequence is an R-minimizing solution. Moreover if the R-minimizing solution uy is unique then uk ! uy with respect to sU . Further we recall some classical results from functional analysis. Lemma 1. The norms in Banach spaces are weakly lower semicontinuous. Theorem 5 (Kakutani). Let X be a Banach space. The following are equivalent: (i) X is reflexive (ii) The closed unit ball of X is weakly compact. (iii) Every bounded set admits a weakly convergent subsequence. Similar statement is due to [17]: Theorem 6. A normed linear space is reflexive if its unit ball is weakly compact. Further we rigorously define the BV space. Details about Hk ðXÞ can be found e.g. in [14]. 4.2. Space of BV functions 1 n n We will frequently make use of functions of bounded variation. We denote by C 1 c ðX; R Þ the set of functions in C ðR Þ with compact support in X.
11588
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Definition 2. By BVðXÞ we denote the subspace of functions u 2 L1 ðXÞ such that the following quantity is finite
Jðu; XÞ ¼ sup
Z X
n uðxÞdivnðxÞdx j n 2 C 1 c ðX; R Þ; knkL1 ðX;Rn Þ 6 1 :
We define the norm kukBVðXÞ ¼ kukL1 ðXÞ þ Jðu; XÞ. We point out that BVðXÞ endowed with the norm kukBVðXÞ is a Banach space and when X R2 we have that BVðXÞ L2 ðXÞ. From [1] we have the following results. Lemma 2. The norm k kBV ðXÞ is weakly lower semi-continuous with respect to the Lp topology for 1 6 p < 1. Lemma 3. Let S be a BV-bounded set of functions. Then S is relatively compact in Lp ðXÞ for 1 6 p < n=ðn 1Þ. S is bounded and thus relatively weakly compact in Lp ðXÞ for p ¼ n=ðn 1Þ and n P 2. 5. Domain decomposition method In the regularization theory, several regularizations are widely used to treat the ill-posedness of the underlying inverse problems. For the clear demonstration of our ideas, let us choose two particular regularizations: the L regularization and the total variation regularization, often called also the bounded variation (BV) regularization. Both regularizations have their drawbacks. The L regularization suffers from the oversmoothing when a discontinuous solution has to be recovered. The BV regularization appears to detect the discontinuities very well, however, when a smooth solution is to be recovered, a staircase effect appears. Therefore, if the unknown solution has discontinuities and on some parts of the domain is smooth, none of the two mentioned regularizations does recover this solution properly. As already mentioned in Introduction we follow a very simple idea: Let us use the L regularization on that part of the domain where we suppose the solution to be smooth and the BV regularization on that part of the domain where the discontinuities are likely to appear. In Section 4 we built up a theoretical framework for incorporation of several regularizations. The application we had on mind was the one of combination of several different regularizations on different parts of the computational domain. We rigorously define our inverse problem. Consider a domain X 2 Rn that is split into two non-overlapping subdomains X1 ; X2 so that X ¼ X1 [ X2 ; m ¼ 2. The first subdomain represents the part, where the solution is expected to be smooth, and the second subdomain should cover the discontinuities of the solution. At this point we do not exploit how this domain decomposition can be achieved. This topic is discussed in Section 6. The forward problem is represented by an (not necessarily linear) operator F operating between functional spaces U and V. For example, for the image recovery problems it is the identity operator, for the inverse problems in diffusion it can be an operator giving a solution of a partial differential equation. The only assumptions we put on F are those from Assumption 1. Further take U 1 ¼ Hk ðX1 Þ and U 2 ¼ BVðX2 Þ, where k is the desired order of expected regularity of the solution. We construct U by
U ¼ fu 2 Ln=ðn1Þ ðXÞ : ujX1 2 U 1 ; ujX2 2 U 2 g:
ð11Þ
Introducing the projector Pi : U ! U i being the restriction on Xi for i ¼ 1; 2 we have that U ¼ U 1 U 2 . We equip U with the following norm kukU ¼ kP1 ukHk ðX1 Þ þ kP2 ukBVðX2 Þ . Further take V ¼ L2 ðX1 Þ Lp ðX2 Þ, where p ¼ n=ðn 1Þ, the norm in V being kv kV ¼ kP 1 v kL2 ðX1 Þ þ kP2 v kLp ðX2 Þ . Consider the regularized least-square solution of the operator equation (6) with the noisy data v d for r ¼ 2
minfkFðuÞ v d k2V þ ab1 kP 1 ukHk ðX1 Þ þ ab2 kP2 ukBVðX2 Þ g: u2U
ð12Þ
Notice, that in the specific case F ¼ I, the minimization problem (12) is equivalent to two independent minimization problems: on X1 with Hk ðX1 Þ-regularization and on X2 with BVðX2 Þ-regularization. To analyse the minimization problem given by (12) we will apply the theoretical framework from Section 4. If the problem fulfills Assumption 1, all the theoretical results from Section 4 apply automatically. We verify Assumption 1 in the following subsection. 5.1. Verification of Assumption 1 (1) Consider U; U 1 ; U 2 as in (11). Let sV be topology generated by the weak norm on the direct sum L2 ðX1 Þ Lp ðX2 Þ. More specifically, the convergence of v k to v in sV means that P1 v k converges weakly to P 1 v in L2 ðX1 Þ and that P 2 v k converges weakly to P2 v in Lp ðX2 Þ. By sU1 denote the topology generated by the weak norm on Hk ðX1 Þ, by sU2 denote
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
the topology generated by the weak norm on Lp ðX2 Þ. So we have U and V the Banach spaces equipped with the topologies sU ; topologies. (2) Assume that v k ! v in V with respect to sV . By the definition it means that
P1 v k * P1 v P2 v k * P2 v
11589
sV which are weaker than the norm
in L2 ðX1 Þ; in Lp ðX2 Þ:
From Lemma 1, the strong norms in L2 ðX1 Þ and Lp ðX2 Þ are weakly lower semi-continuous. That results in
kP1 v kL2 ðX1 Þ 6 lim inf kP 1 v k kL2 ðX1 Þ ; kP2 v kLp ðX2 Þ 6 lim inf kP2 v k kLp ðX2 Þ : Putting it together we conclude
kv kV ¼ kP1 v kL2 ðX1 Þ þ kP2 v kLp ðX2 Þ 6 lim inf kP1 v k kL2 ðX1 Þ þ lim inf kP2 v k kLp ðX2 Þ 6 lim inf kv k kV ; which verifies that k kV is sequentially lower semi-continuous with respect to sV . (3) This assumption will be verified here only for the case of the identity operator F ¼ I with DðFÞ ¼ Hk ðX1 Þ BVðX2 Þ. We have to verify the continuity of the identity from DðFÞ to V ¼ L2 ðX1 Þ Lp ðX2 Þ with respect to weak topologies. Suppose uk ! u in sU . That means that
P1 uk * P 1 u in Hk ðX1 Þ; P2 uk * P 2 u in Lp ðX2 Þ: It directly follows that also
P1 uk * P1 u in L2 ðX1 Þ: This gives the continuity of the identity operator from DðFÞ to V with respect to the topologies sU and sV . If F is not the identity operator, we have to verify its continuity for each specific case separately. (4) Setting R1 ðuÞ ¼ kukHk ðX1 Þ and R2 ðuÞ ¼ kukBVðX2 Þ we directly have that RðuÞ ¼ b1 R1 ðP1 uÞ þ b2 R2 ðP2 uÞ is proper and convex. From Lemma 1 we know that R1 is sU1 lower semi-continuous. Lower semi-continuity of R2 with respect to sU 2 follows from Lemma 2. (5) Since DðFÞ is a direct sum of closed spaces, DðFÞ is closed too. (6) Suppose fuk g 2 Ma ðMÞ which means
kFðuk ÞkV þ kP 1 uk kHk ðX1 Þ þ kP2 uk kBVðX2 Þ 6 C: For reflexive space Hk ðX1 Þ we know from Theorem 5 that the unit ball is weakly compact. Since fP 1 uk g is bounded in Hk ðX1 Þ then there exists weakly convergent subsequence again denoted by fP 1 uk g. So we have that there exists u1 2 Hk ðX1 Þ such that P1 uk ! u1 with respect to sU1 topology. From Lemma 3 and the boundedness of fP 2 uk g in BVðX2 Þ we know that there exists subsequence again denoted by P2 uk convergent in sU2 topology with the limit u2 2 BVðX2 Þ. We construct u 2 U such that P 1 u ¼ u1 and P 2 u ¼ u2 and as a conclusion of the above preliminaries we state that from the sequence fuk g we are able to extract a subsequence convergent to u in sU topology. We successfully verified Assumption 1 for our setting. The parts (1), (2), (4)–(6) of Assumption 1 have been verified independently of F. The part (3) strongly depends on F and we verified it for the case of the identity operator F ¼ I. So we are able to formulate well-posedness and convergence results for general nonlinear F in the following theorem. Theorem 7. Consider the minimization problem (12). If F : DðFÞ U ! V is continuous with respect to the topologies sU and sV then Theorems 1–4 are valid for this minimization problem. Remark 3. It is possible to consider the BV semi-norm Jðu; X2 Þ instead of the BV norm kukBVðX2 Þ in all our considerations. The statements of Theorems 1–4 will remain valid after one more assumption on F, namely that F does not annihilate constant functions on X2 . The verification of the statements is analogous and will copy the reasoning in [1]. 6. Implementation In this section, we apply the ideas and the results obtained in Sections 4 and 5 to image denoising. The theoretical considerations from Section 5 are starting from the assumption that the image domain X is already split into two regions X1 ; X2 where X1 contains the smooth parts of the image and X2 the discontinuities. To determine the locations of the edges in the image, one can run any edge-detection algorithm that localizes steep variations of the image gray level. For an overview we
11590
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
refer to [20] and the references therein. Any of these techniques include an elaborated mechanism to choose the threshold according to which the decision is made whether edges are present or not at an image point. The lower the threshold, the more edges will be detected, and the result will be increasingly susceptible to noise, and also to picking out irrelevant features from the image. Conversely a high threshold may miss subtle edges, or result in fragmented edges. The aim of this paper is not to design an effective algorithm for edge detection but rather to design an algorithm that can denoise the image after we have the edges detected. Therefore we have focused our theoretical consideration on the denoising part having the image decomposed. For the actual decomposition we use the first order technique applied on the BV regularized image. We choose the BV-regularized image for the actual detection of the edges since the edges remain preserved after the BV regularization. Unlike in Sections 4 and 5 we use the BV seminorms Jðu; X2 Þ; Jðu; XÞ instead of the BV norms kukBVðX2 Þ ; kukBVðXÞ as justified in Remark 3. The forward operator is F ¼ I, the original image is v. The image v is polluted by noise of a known level d. The polluted image is denoted by v d . The ultimate goal of image denoising is to recover the original image v from v d . Obviously this is not possible exactly. We try to approximate v by a minimizer of the regularized cost functional (7). The quality of the reconstruction is determined by two factors: a and the ratio b1 =b2 . The parameter a in (7) determines the global weight of the regularization term RðuÞ. On the other hand, the parameters b1 ; b2 in (8) determine the relative weight of particular regularizations, in our case the Hk ðX1 Þ and the BVðX2 Þ regularizations. Consequently, the Hk regularization is weighted by c1 :¼ ab1 and the BV regularization is weighted by c2 :¼ ab2 . Given that c1 ; c2 are determined, Theorem 4 yields the convergence of a general method, satisfying the assumptions of the theorem, to the best-approximate solution in the sense of Definition 1. We minimize all the involved cost functionals by conjugate gradient method. We propose the following algorithm. The first and the second step comprise the decomposition of the domain and the actual denoising is done in the third step. Algorithm 1 1. Find u 2 BVðXÞ as a minimizer of
T BV ðuÞ :¼ kFðuÞ v d k2V þ cBV Jðu; XÞ:
ð13Þ
The regularization constant cBV will be determined by the discrepancy principle [16], see below. 2. Find the edges of u and decompose the image domain to X2 which contains all the edges and X1 :¼ X n X2 . 3. Find w 2 U as a minimizer of
T MIXED ðwÞ :¼ kFðwÞ v d k2V þ cMIXED kwkHk ðX1 Þ þ cBV Jðu; X2 Þ;
ð14Þ
where cMIXED is again determined via the discrepancy principle. 6.1. Discrepancy principle The discrepancy principle is a classical criterion used to determine regularization constant for ill-posed problems [9,16]. Solving FðuÞ ¼ v d is feasible only up to the level of noise, i.e we should not expect to find a regularized solution uda such that kFðuda Þ v d k < d (a is the corresponding regularization constant). In the same time one wants to introduce as much regularization as possible. When applying the discrepancy principle a is successively decreased untill kFðuda Þ v d k < Dd. The constant D 1 determines how much regularization is introduced and strongly influences the result. We however found the successive increase of a to be numerically more stable. Putting this together we propose the following simple update rule for a. 1. Set a to a small constant. 2. While kFðudak Þ v d k < Dd set akþ1 ¼ 2ak . We can get a more accurately by bisectioning on interval ak1 ; ak ½ for the final index k. 6.2. Domain decomposition To obtain the decomposition of the domain X into X1 and X2 from the step 2 of Algorithm 1 we employ the following simple approach. We compute the norm of the gradient of u. If it exceeds a certain value for a point of the image, this point belongs to X2 . Except this point we take into X2 also the points from some immediate neighborhood. The rest of image points belongs to X1 . 7. Numerical experiments We apply Algorithm 1 to eight real world 256256 images depicted in Fig. 1. Real images usually contain flat regions, smooth parts and textures. Usually none of these components vanishes at any point of an image. However, we assumed that we had an a priori information that some components are vanishing on a part of the domain. Namely we assume that the
11591
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Fig. 1. Images.
Table 1 Regularization constants and relative errors of L, of BV and of MIXED regularizations. Image
Caesar
Clown
Deer
Girl
Collection
Lighthouse
Lena
Airport
eL eBV eMIXED
0.0923 0.1128 0.0842 0.84 9.6 1.45
0.1047 0.1199 0.0876 0.675 8.4 2.657
0.1667 0.1411 0.1669 0.2225 1.63 0.888
0.0902 0.1094 0.0701 1.05 17.6 2.3133
0.1057 0.0904 0.1006 0.4275 3.22 0.7882
0.1009 0.1336 0.0952 1.6 29.28 2.87
0.1182 0.1272 0.0996 0.3975 3.8 3.7266
0.1085 0.1387 0.0997 1.39 24.96 6.6459
cL cBV cMIXED
images are composed of smooth regions separated from each other by sharp edges. This does not hold for general images. Consequently, one cannot expect that the domain decomposition of Algorithm 1 will be more powerful for general images than some algorithms from the Introduction. However, we will demonstrate that the domain decomposition works well. The discrepancy principle stopping constant D is set manually to DBV ¼ 0:5 for the BV regularization and DL ¼ DMIXED ¼ 1:3 for the smooth and the mixed regularizations (the regularization parameter cL of the smooth regularization is obtained in the same way as cBV by discrepancy principle). We use these values for all the experiments. First, we apply the uniform noise of 10% to all the images. In Table 1 one can find the corresponding numerical results: eL , eBV , and eMIXED are the relative L2 errors between the original image and the denoised one for the L, the BV and the MIXED regularization, respectively. cL ; cBV and cMIXED are the corresponding regularization parameters obtained by Algorithm 1. In majority of the cases the mixed regularization performs better than other regularizations and achieves the smallest relative error with respect to the original image. It is however misleading to compare the results on the basis of the relative errors only. Even if the reconstruction can be very close in L2 -sense, small details can be modified in a way that the image is visually degraded. For all the images, Algorithm 1 finds bigger cMIXED than cL even if DL ¼ DMIXED . It is caused by the following: If the L regularization is applied everywhere in the domain, it smears the edges. If cL was bigger than it is, the edges would be smeared much more and the reconstruction would be less accurate. Let us recall that discrepancy principle finds the maximal c’s such that the reconstructions are as close as possible to the noisy image, i.e. as one expects also as close as possible to the original one - up to the level of noise. When the L regularization is applied only on the smooth part of the domain X1 , it does not smear the edges in X2 . More L regularization can be applied in the smooth areas of the image and the reconstruction will be still close to the noisy image. The noise is then removed much better. Let us discuss three different images more closely. In Fig. 2 we present the L, the BV and the MIXED reconstructions for the Caesar image. The MIXED reconstruction is substantially better than the rest. First, the edges are recovered sharper than by the BV regularization (notice that the regularization constant cBV is the same). In the case of the MIXED regularization the BV
11592
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Fig. 2. Caesar: (a) original; (b) L-reconstruction; (c) BV-reconstruction; (d) MIXED-reconstruction.
Fig. 3. Caesar: (a) mask; (b) residuum of L-reconstruction; (c) residuum of BV-reconstruction; (d) residuum of MIXED-reconstruction.
regularization is applied only in some immediate neighborhood of the edges. Here the image is closely approximated by a function with a small BV norm - two flat regions divided by the edge. The gray scales are also recovered more accurately. The points in the smooth part of the domain X1 do not influence the intensities of points in X2 . The staircase effect is not visible. On the other hand, as already explained, more L regularization is applied in X1 . The noise is much better removed then in the case of plain L regularization. Fig. 3(a) displays the decomposition of X into X1 and X2 . White color represents the edges. i.e. X2 . Fig. 3(b)–(d) display the residua of the L, BV and MIXED reconstructions. The residua are normalized for clarity. The MIXED residuum contains clearly the least information, notice the details around Caesar’s nose and mouth.
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
11593
Fig. 4. Airport: (a) original; (b) L-reconstruction; (c) BV-reconstruction; (d) MIXED-reconstruction.
Fig. 5. Airport: (a) mask; (b) noisy image.
In Fig. 4 the results for the airport image are presented. This image contains a lot of small details. The MIXED reconstruction is clearly the best one. Even small details are recovered very well. Edges are sharp and intensities are correct. The corresponding decomposition of X and the noisy image are presented in Fig. 5.
11594
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Fig. 6. Lena: (a) original; (b) L-reconstruction; (c) BV-reconstruction; (d) MIXED-reconstruction.
Fig. 7. Lena: (a) mask; (b) residuum of L-reconstruction; (c) residuum of BV-reconstruction; (d) residuum of MIXED-reconstruction.
At last we present the results for Lena image. From Fig. 6 we can again conclude that the MIXED regularization performs best. The result seems a bit oversmoothed, e.g. the details of the feathers and on the face are lost. However the noise is removed much better then in the case of L regularization and the edges are sharper then for plain BV regularization. The corresponding domain decomposition and the residua are depicted in Fig. 7. The recovered images of Lena by MIXED regularization for different noise levels are depicted in Fig. 8. The respective numerical results are to find in Table 2.
11595
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Fig. 8. Lena: results obtained by mixed regularization for different noise levels of (a) 1%; (b) 3%; (c) 5%; (d) 20%.
Table 2 Lena: results obtained by mixed regularization for different noise levels of (a) 1%; (b) 3%; (c) 5%; (d) 20%. noise
1%
3%
5%
20%
eL eBV eMIXED
0.0158 0.0178 0.0138 0.0472 0.54 0.7786
0.0422 0.0481 0.0355 0.1838 2.06 9.1875
0.0635 0.0751 0.0541 0.3875 4.92 20.5375
0.1987 0.2642 0.1964 10.56 135.68 15.1287
cL cBV cMIXED
8. Conclusions We presented a theoretical framework for mixed Tikhonov regularization in Banach spaces. We first extended the results of [12], where the Tikhonov regularization in Banach spaces was tackled through completeness in a weaker than the norm topology. Particularly we looked for a solution in a direct sum of several Banach spaces and the weak topology was defined using projectors on these Banach spaces. Such projectors are naturally defined, if one considers a decomposition of the domain to certain subdomains and the Banach spaces are defined on these subdomains. The results obtained for domain decomposition based Tikhonov regularization are applicable for general possibly non-linear inverse problems, if the operator corresponding to the direct problem is continuous with respect to the weak topologies, see Theorem 7. We applied the framework to image denoising problem, when the direct problem is simply described by identity map. The image was decomposed to two parts. On one part, the smoothing L regularization was applied, while on the second part the BV regularization. We showed that the image denoising problem fulfills all the assumptions of the theoretical framework.
11596
I. Cimrák, V. Melicher / Applied Mathematics and Computation 218 (2012) 11583–11596
Finally, we implemented an algorithm which performs the mixed L–BV regularization. The results for a set of real world images confirm that the domain decomposition idea works well. Acknowledgment Large part of the presented work was done during the period when Ivan Cimrák was supported by the Fund for Scientific Research – Flanders FWO, Belgium. References [1] R. Acar, C.R. Vogel, Analysis of bounded variation penalty methods for ill-posed problems, Inverse Problems 10 (1994) 1217–1229. [2] G. Aubert, L. Vese, A variational method in image recovery, SIAM Journal of Numerical Analysis 34 (1997) 1948–1979. [3] J. Aujol, G. Gilboa, T. Chan, S. Osher, Structure-texture image decomposition–modeling, algorithms, and parameter selection, International Journal of Computer Vision 67 (2006) 111–136. [4] J.-F. Aujol, G. Aubert, L. Blanc-Féraud, A. Chambolle, Image decomposition into a bounded variation component and an oscillating component, Journal of Mathematical Imaging and Vision 22 (2005) 71–88. [5] P. Blomgren, T.F. Chan, P. Mulet, Extensions to total variation denoising, in: F.T. Luk (Ed.), Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, vol. 3162, 1997, pp. 367–375. [6] T. Chan, A. Marquina, P. Mulet, High-order total variation-based image restoration, SIAM Journal on Scientific Computing 22 (2000) 503–516. electronic. [7] P. Charbonnier, Reconstruction d’image: regularisation avec prise en compte des discontinuitis, Ph.D Thesis, University of Nice-Sophia Antipolis, 1994. [8] I. Cimrák, Analysis of the bounded variation and the G regularization for nonlinear inverse problems, Mathematical Methods in the Applied Sciences 33 (2010) 1102–1111. [9] H.W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Mathematics and its Applications, vol. 375, Kluwer Academic Publishers, Dordrecht, 1996. [10] G. Gilboa, N. Sochen, Y. Zeevi, Variational denoising of partly-textured images by spatially varying constraints, IEEE Transactions on Image Processing 15 (2006) 2281–2289. [11] J. Gilles, Noisy image decomposition: a new structure texture and noise model based on local adaptivity, Journal of Mathematical Imaging and Vision 28 (2007) 285–295. [12] B. Hofmann, B. Kaltenbacher, C. Pöschl, O. Scherzer, A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators, Inverse Problems 23 (2007) 987–1010. [13] K. Ito, K. Kunisch, Bv-type regularization methods for convoluted objects with edge flat and grey scales, Inverse Problems 16 (2000) 909–928. [14] A. Kufner, O. John, S. Fucˇík, Function spaces, Monographs and Textbooks on Mechanics of Solids and Fluids, Mechanics: Analysis, Noordhoff International Publishing, Leyden, 1977. Academia, Prague. [15] Y. Meyer, Oscillating patterns in image processing and in some nonlinear evolution equations, March 2001. The Fifteenth Dean Jacquelines B. Lewis Memorial Lectures. [16] V.A. Morozov, On the solution of functional equations by the method of regularization, Doklady Akademii Nauk SSSR 7 (1966) 414–417 (in Russian). [17] D.M. Oberlin, A measure-theoretic proof of a theorem on reflexivity, Proceedings of the American Mathematical Society 41 (1973) 325–326. [18] L. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1992) 259–268. [19] D.M. Strong, P. Blomgren, T.F. Chan, Spatially adaptive local feature-driven total variation minimizing image restoration, in: F. Preteux, J.L. Davidson, E.R. Dougherty (Eds.), Statistical and Stochastic Methods in Image Processing II, vol. 3167 of Proceedings of Society of Photo-Optical Instrumentation Engineers (SPIE), Bellingham, WA, 1997, Soc Photo Opt Instrumentation Engineers, SPIE – INT SOC OPTICAL ENGINEERING, pp. 222–233, Conference on Statistical and Stochastic Methods in Image Processing II, San Diego, CA, July 31– August 01, 1997. [20] D. Ziou, S. Tabbone, Edge detection techniques – an overview, International Journal of Pattern Recognition and Image Analysis 8 (1998) 537–559.
