RETINEX BASED ON EXPONENT-TYPE TOTAL VARIATION SCHEME

. Retinex theory deals with compensation for illumination eﬀects in images, which has a number of applications including Retinex illusion, med- ical image intensity inhomogeneity and color image shadow eﬀect etc.. Such ill-posed problem has been studied by researchers for decades. However, most exiting methods paid little attention to the noises contained in the images and lost eﬀectiveness when the noises increase. The main aim of this paper is to present a general Retinex model to eﬀectively and robustly restore images degenerated by both illusion and noises. We propose a novel variational model by incorporating appropriate regularization technique for the reﬂectance component and illumination component accordingly. Although the proposed model is non-convex, we prove the existence of the minimizers theoretically. Furthermore, we design a fast and eﬃcient alternating minimization algorithm for the proposed model, where all subproblems have the closed-form solutions. Applications of the algorithm to various gray images and color images with noises of diﬀerent distributions yield promising results. minimization algorithm.


1.
Introduction. Decomposition an image into meaningful components is the real substance of many tasks in image processing, such as image denoising, which decomposes the image into signal parts and noise parts [32,7]; image structure-texture decomposition, where images are modeled as a combination of geometrical information and textural information [1,2]. Retinex [19] is a related but different topic, which was originally proposed by Land and McCann [19] as a model of color perception of the human visual system. It reveals that human vision tends to see the same color in a given image regardless the light and the color of objects remains relatively constant under varying illumination. Thus, the goal of Retinex is to decompose the illumination from the reflectance in given images. In fact, such problem also exists in medical image processing, the so-called intensity inhomogeneity, which leads to the intensity variation even for the same tissue within an image. The importance of Retinex is that it can facilitate other image analysis techniques such as segmentation, registration, etc., that relying on the assumption of uniform intensity.
To set up the Retinex problem mathematically, we focus on decomposing a given image I into the reflectance component R and the illumination component L as The following general assumptions exists as the basic laws for Retinex: • The reflectance R(x, y) is a piecewise constant function; • The illumination L(x, y) is a spatial smooth function. Since the logarithmic transform can preserve image edges and simplify the relationship between R and L, most Retinex algorithms are proposed in the logarithmic domain, which gives i(x, y) = r(x, y) + l(x, y), (2) where i = log(I), l = log(L), r = log(R).
Various implementations and algorithms have been studied for Retinex problem. Path-based methods were originated by Land and McCann [19,20], and further studied in [10,31]. Recursive algorithms were proposed by using a recursive matrix calculation to replace the path computation in [21,26]. Homomorphic filtering type Retinex algorithms modeled the reflectance as a low-pass version of the given image based on a convolution with a wide Gaussian kernel [14,17]. Bertalmío et al. [4] proposed a kernel-based Retinex algorithm relying on the computation of the expectation value of a suitable random variable weighted with a kernel function. Retinex was also considered in an alternative perspective [5,29], the so-called color correction, both of which were built up in the variational framework.
In PDE-based models and variational formulations, Morel et al. [27] formulated the illumination as a spatially smooth image and the reflectance as a piecewise constant image, which was solved as a Poisson equation by Fast Fourier Transform (FFT). Based on the assumption of the spatial smoothness of the illumination, Kimmel et al. [18] proposed a variational model by penalizing on the reflectance. Elad [13] proposed a noniterative Retinex algorithm by utilizing two bilateral filters as the regularization term. Ng and Wang [28] proposed a similar model by penalizing both the reflectance and illumination images. Ma and Osher [25] developed a Total Variation (TV) [32] based model to extract the reflectance image with a data term in gradient field. Zosso et al. [37] further extended the TV-based Retinex models to a unified non-local formulation. Chang et al. [8] proposed a variational model for Retinex by representing the reflectance component with more details by a learned dictionary. Duan et al. [11,12] investigated variational models for Retinex using 0 quasi-norm to regularize the reflectance image. Liang and Zhang [23] proposed a convex Retinex model, which decompose the gradient field of images into salient edges and relatively smoother illumination field. Various algorithms have gained great success in dealing with Retinex problems. However, these methods paid little attention to the noises contained in given images, as they are all based on the image model (1), which does not model noises as a component. Especially, when the logarithmic transform is taken as (2), the distribution of the noises becomes more complicated, and even harder to remove. In fact, illusion (intensity inhomogeneity) and noises may simultaneously exist in given images, such as MR images, ultrasound images etc.. Therefore, when the observed images are corrupted by noises, these methods lost their effectiveness in recovering the reflectance image with high quality.
1.1. Our contributions. In this paper, we propose a general variational model for images in presence of both illusion and noises, that is we consider the following image model (3) I(x, y) = (R(x, y) · L(x, y)) ⊕ n(x, y).
The multiplicative relation between R and L can be transformed into the addition one as follows where v = log V , ⊕ represents the effect of noises, which indicates + for additive noise and × for multiplicative noise. Based on (4), we reported in [24] the idea of recovering the noise-free image V from the given image I, and decomposing it into the reflectance r and illumination l in the logarithmic domain. However, only preliminary results were presented and some advantages were observed, where we extend the idea in both theoretically and numerically in this work. To be specific, we use total variation as the regularization for V and choose the date fidelity according to the distribution of the noise model. For the decomposition, the first-and secondorder TV regularization are implemented as [23]. Furthermore, we propose an efficient alternating minimization (AM) algorithm, all subproblems of which can be efficiently solved by the closed-form solutions. Numerous experiments on synthetic and real images demonstrate that our model gives better results in comparison with other state-of-the-art Retinex methods, especially when noises is significant. In addition, three main contributions are presented: 1. We propose a variational model to recover images corrupted by illusion and noises, which is guaranteed with the existence of minimizers. 2. Our proposed model is capable to deal with different types of noises, such as Gaussian noise, impulsive noise [6] and Poisson noise [22] etc., as long as we modify the data fidelity term according to the noise distribution. 3. Our framework can be easily extended to other Retinex models, such as Ng and Wang's Retinex model [28], 0 regularized Mumford-Shah model [11,12], etc.. The rest of the paper is organized as follows. In Section 2, we introduce the notations that will be used in the rest of the paper. We briefly review the variational models for Retinex. In Section 3, we propose our new model for recovering images corrupted by illusion and noises. Furthermore, the existence of the solution to the proposed model is provided and the proofs of the mathematical results are shown. In section 4, we design an alternating minimization algorithm for the proposed model and discuss the solutions to each sub-minimization problem. In Section 5, we present the numerical experiments on both synthetic medical images, visual illusion images and real medical images, color images. We then conclude the paper in Section 6 with some final remarks and future prospects.
2. Setting and previous works.

2.1.
Discretization. Without loss of generality, we present a grayscale image as a two dimension matrix of size N × N , the size of matrix here can also be M × N . We denote X the Euclidean space R N ×N , which is equipped with the usual inner product and Euclidean norm as ·, · and · 2 , respectively. The discrete gradient operator is a mapping ∇ : X → Y, where Y = X × X. For u ∈ X, ∇u is given by where j, k = 1, . . . , N . Here we useD + x andD + y to denote forward difference operators with the periodic boundary conditions.
We also equip the space Y with the inner product p, q = p 1 , q 1 + p 2 , q 2 and the norm p 1 = 1≤i,j≤N whereD − x andD − y to denote backward difference operators with the periodic boundary conditions. 2.2. TVH1 model. We review some important variational models for Retinex. Kimmel et al. [18] proposed a variational Retinex model as follows where α and β are positive constants. The first term in (5) is used for the fidelity, while the reflectance part and the illumination part are regularized by the TV norm and H 1 norm, respectively.
Ng and Wang further studied the total variation model for Retinex in [28], where they explicitly minimized both r and l in the energy functional with box constraints min r≤0,l≥i where τ is a positive constant, and the constraint l ≥ i is given based on the assumption that R is the reflectivity and 0 < R ≤ 1. [23] presented a high order TV based model for reflectance and illumination decomposition, which reads min r∈Br,l∈B l

HoTVL1 model. Recently, Liang and Zhang
where the first-and second-order TV are introduced to regularize the reflectance and illumination, respectively. The box constraints are imposed on both r and l, which vary with different applications.
2.4. L0MS model. Duan et al. [12] built up an unconstrained variational model based on the Mumford-Shah formulation for images with intensity inhomogeneity as follows , and for a given pixel p ∈ Ω, N p = {q| |p − q| ≤ ρ} is its neighborhood with a radius ρ. In (8), r and l are regularized by the 0 quasi-norm and the H 1 norm, respectively. Unlike other methods, L0MS model gives an exact piecewise constant solution of r, which is a preliminary results for structural segmentation. Thus, the 0 regularized Retinex model has been reformulated for three-dimensional applications based on the highorder regularization in [9].
3.1. The proposed model. In real applications, the observed images are usually contaminated by both illusion and noises. Consequently, we aim to recover the noise-free image, and decompose it into reflectance image and illumination image. More specifically, we propose a novel Retinex model by combining the HoTVL1 model (7) and a data fidelity term deriving from the noise models. We use the additive Gaussian noise as an example, and propose the following unconstrained minimization problem where V = e v according to the image model (4). The main advantage of using the exponential function other than the logarithmic function is its convexity [16]. Obviously, the above model can be easily adapted to other noise models by choosing an appropriate data fidelity term, which will be discussed in the section of numerical experiments.

Existence of solution.
Although the energy functional F(v, r, l) in (9) is non-convex, we are still able to show the existence of its minimizers. On the first place, we prove that the objective function F(v, r, l) is coercive as follows.
Proof. The proof is motivated by Lemma 3.8 in [16]. The lower bound of the discrete TV is given by where L is the first-order difference matrix defined by the one-sided difference matrix on the horizontal direction D x and the vertical direction D y as The null space of L is the set {c 1 1} with c 1 being a scalar. On the other hand, the lower bound of the discrete high-order TV is given as where H is defined as with D xx , D xy , D yx and D yy being the second-order difference matrixes. The null space of H is the set {c 2 1} with c 2 being a scalar.
We consider two cases: For case (i), consider (v, r, l) / ∈ Γ with (v, r, l) 2 → ∞, we can discuss according to the value of v, r and l and easily obtain For case (ii), consider (v, r, l) ∈ Γ, there is v = r + l, r = c 1 1, and l = c 2 1. Since (v, r, l) 2 → ∞, i.e., at least |c 1 | → ∞ or |c 2 | → ∞,we have In order to establish the existence of the minimizers for (9), we first present some preliminaries for the bounded Hessian (BH) space. Let Ω be an open subset of R N ×N with Lipschitz boundary. Then we set BV(Ω) = {u| |u| T V < ∞}, where |u| T V = |Du|. We consider v, r in the bounded variation space BV(Ω) and l in the bounded Hessian space BH(Ω), where BH(Ω) = {l| D 2 l 1 < ∞}. We extend some main properties of the high-order total variation space in [23,30] to the following forms of the discretized space BH(Ω).
• Suppose that {u k } k∈N is bounded in BH(Ω), then there exists a subsequence {u kj } j∈N and u ∈ BH(Ω) such that {u kj } j∈N weakly * converges to u.
• If Ω has a Lipschitz boundary and it is connected, then it can be shown that there exists positive constants C 1 , C 2 such that Theorem 3.2. The problem (9) has at least one solution (v * , r * , l * ) ∈ BV(Ω) × BV(Ω) × BH(Ω).
Proof. Let us pick a minimizing sequence (v n , r n , up to a subsequence also denoted by {r n } after relabeling, there exists a function r * ∈ BV(Ω) such that (a) r n → r * strongly in 1 (Ω), (b) r n → r * a.e. x ∈ Ω, (c) ∇r n ∇r * in the sense of measure. The lower semi-continuity of the total variation leads to (11) ∇r * 1 ≤ lim inf n→∞ ∇r n 1 .
Similarly, we can derive that {l n } is bounded in 1 (Ω). By virtue of the embedding inequality (10), we also have ∇l n 1 < C 1 D 2 l n 1 + C 2 l n 1 < M , where M is a constant number for all n ∈ N. Thus {l n } is uniformly bounded in BH(Ω), up to a subsequence also denoted by {l n } after relabeling, there exists a function l * ∈ BH(Ω) such that (a) l n → l * strongly in 1 (Ω), (b) l n → l * a.e.

BV(Ω)
is uniformly bounded. So e v n BV(Ω) is uniformly bounded. The proof of v is the same as r, up to a subsequence also denoted by {v n } after relabeling, there exists a function v * ∈ BV(Ω) such that v n → v * a.e. x ∈ Ω. Then, by the lower semi-continuity of the total variation and 2 -norm, we also have Finally, since r n → r * a.e. x ∈ Ω, l n → l * a.e. x ∈ Ω, v n → v * a.e. x ∈ Ω, Fatou's lemma gives that Combining inequalities from (11) to (14), we deduce that which indicates that (v * , r * , l * ) is a minimizer of problem (9). This completes the proof.
4. The alternating minimization algorithm. In this section, we discuss the minimization of (9). Since the variables v, r and l are coupled together, an accurate joint minimization can be costly. Thus, we use the alternating minimization algorithm to solve the minimization problem (9), through which one can obtain the minimizer approximately [33,15,34]; see Algorithm I.

4.1.
Sub-minimization problem with respect to v. In (15), the gradient and exponential function are coupled together, which increases the difficulty of computation. Therefore, we introduce a new variable u = e v , and reformulate the minimization problem (15) as the following constrained minimization problem min v,u We are going to minimize the variable u, which belongs to the category of 1regularized problems with the general form as follows where both Φ(u) 1 and H(u) are convex functions. The key idea to efficiently solve the minimization problem (19) is to split the two portions in (19) as [15,34]. For our case, rather than considering (18), we introduce another auxiliary variable p and consider the following constrained optimization problem min v,u,p Based on the augmented Lagrangian method, we can obtain the equivalent unconstrained optimization problem as where Λ 1 and Λ 2 are Lagrange multipliers, and r 1 , r 2 are positive penalty parameters. It is straightforward to separate (21) into the following three subproblems The subproblems (22), (23) and (24) can be efficiently solved based on closedform solutions. For (22), its Euler-Lagrange equation gives a nonlinear equation due to the exponential function of v. Let g(v) = Λ n 1 , −e v + r1 2 u n − e v 2 2 , which is linearized by its Taylor expansion as follows [3] By adding an extra proximal term [36] and deleting the constant terms, we arrive at the following minimization problem for v the solution of which gives us v n+1 = Λ n 1 e v n − r 1 e v n (e v n − u n ) + v n + µ(r n + l n ) 1 + µ .
For (23), its optimality condition gives us the following linear equation where I is the identity operator. By the periodic boundary condition for u, (26) becomes a block circulant system, which can be efficiently solved by the fast Fourier transform (FFT). Denoting F(u) as its Fourier transform, we can express the solution to u as follows: where p n = (p n,1 , p n,2 ), Λ n 2 = (Λ n,1 2 , Λ n,2 2 ). For (24), we can efficiently solve it using the shrinkage operator as follows where shrink(ω, λ) = ω ω 1 * max( ω 1 − λ, 0). Besides, we update the Lagrange multipliers Λ 1 and Λ 2 according to their gradient ascent direction as Λ n+1

4.2.
Sub-minimization problems with respect to r and l. As observed, both the sub-minimization problems with respect to r and l are 1 regularized optimization problems. Therefore, we can discuss the solutions to (16) and (17) together. Similar to the subproblem of u, we separate the 1 term and 2 term by introducing new variables x = ∇r and y = ∇ 2 l, and solve the corresponding constrained optimization problems based on the augmented Lagrangian method. Here, we omit the implementation details and provide the updating scheme to all variables as follows where Λ 3 and Λ 4 are Lagrange multipliers, and r 3 , r 4 are penalty parameters.
For the subproblem (29) and (31), we can solve them using the FFT under the periodic boundary conditions. For the subproblem (30) and (32), there are the closed-form solutions based on the shrinkage operators. In addition, we update the Lagrange multipliers as follows

Numerical experiments.
In this section, we conduct a series of experiments on synthetic images with different illusions and noises of different distributions to demonstrate the performance of the proposed method.

Tests of gaussian noises.
There are several parameters in our model (9), i.e., the regularization parameters α, β, γ, the penalty parameter µ and the theoretical parameter τ . The regularization parameter α and β are most critical parameters in our model. In the numerical tests, we choose α from α = 0.001 to α = 0.1 according to the severity of noises while choose β from β = 0.001 to β = 0.05 depending on the illusion. For different test images, the parameters γ, µ, τ are fixed as γ = 0.08, µ = 0.7 and τ = 1e − 5. All other penalty parameters are set to r i = 0.2, for i = 1, 2, 3, 4. Three other variational methods, including TVH1 model [28], HoTVL1 model [23] and L0MS model [12], are evaluated and compared with our model. We list the implementation details of each algorithm as follows x TVH1: the model is sensitive to the parameter of the regularization term on reflectance. We adjust β from β = 0.1 to β = 25, and fix α = 0.1, τ = 1e − 5 for different examples. y HoTVL1: the model is sensitive to the parameter of the regularization term on reflectance. We tune α from α = 0.01 to α = 0.1, and set β = 10, τ = 1e − 5.
5.1.1. Stable performance for noises. In magnetic resonance imaging (MRI), the reconstructed images are often corrupted by the bias field due to the inhomogeneous illuminations. We use two slices of T 1 -weighted brain MR image as examples, which are contaminated by 3%, 5%, 7%, 9% noise and 40% intensity nonuniformity (downloaded from

5.1.2.
Stable performance for intensity inhomogeneities. We select one slice from the T 1 -weighted brain volume and generate 10 images by adding Gaussian white noise of mean 0 and variance 0.001 and intensity inhomogeneities of different profiles. We plot the PSNR and MSSIM obtained from the four methods in FIGURE 2, which illustrates that our model gives the best denoising results. Indeed, our model is shown stable with respect to different profiles of intensity inhomogeneities. Moreover, we present two of the 10 examples and the corresponding results in   On the other hand, we evaluate the performance of bias correction using the coefficient of variations (CV), which is defined for each tissue T as  We use the first image in FIGURE 3 as an example, and plot the relative errors in FIGURE 6 (a) and (b). Besides, we also track the numerical energy of the objective function in FIGURE 6 (c). It is shown that both the relative errors of r, l and the numerical energy decay as the iteration increase, which demonstrate that our model converges well numerically. or erroneous transmission [6]. Based on the Bayesian statistic, one can derive a 1 fidelity term, which can exactly fit uncorrupted pixels and perfectly regularize the corrupted pixels by the impulsive noise. Indeed, the 1 fidelity term only affects the sub-minimization problem w.r.t. v, which has the fast numerical solver [35]. In particular, we consider the following minimization problem for images corrupted by salt-and-pepper noise: min v,r,l I − e v 1 + α ∇e v 1 + β ∇r 1 + γ ∇ 2 l 1 + We use the Adelson's checkerboard shadow image and the Logvinenko's cube shadow image as examples, which are two typical test images for Retinex illusion. We add the salt-and-pepper noise from 20% to 50% to both images. The reflectance component R and illumination component L obtained by the proposed method are displayed in FIGURE 9 and FIGURE 10. The visual results demonstrate that our model can well decompose the reflectance and illumination although the test images contain impulsive noises.
For the impulsive noises, we fix β = 0.05, γ = 0.08, µ = 0.6, τ = 1e − 5 and r i = 0.5 for i = 1, · · · , 5, while we select α according to the level of the noises from α = 0.6 to α = 0.9. 5.4. Extension to Poisson noise. The Poisson noise is another common seen noise, which may be contained in radiography, fluorescence microscopy and positronemission-tomography images [22]. According to the the characteristic of Poisson distribution, the Kullback-Leibler (KL) divergence is used as the fidelity term. Therefore, we propose the following minimization problem for images containing Poisson noise: In FIGURE 11, we add Poisson noise into two synthetic images. Both R and L are presented in FIGURE 11,which show that the proposed model can obtain visually preferable results compared to the original ones. For this experiment, we set α = 0.03, β = 0.015, γ = 0.08, µ = 0.6, τ = 1e − 5 and r i = 0.5 for i = 1, · · · , 5.
6. Conclusion and future works. In this paper, we have presented an efficient variational model for Retinex, which was developed based on a new image model by decomposing a given image into reflectance, illumination and noises. The proposed model was shown coercive, which guarantee the existence of the minimizers. We designed an efficient alternating minimization algorithm, where all subproblems can be solved by the closed-form solutions. The framework is easy to be adopted to other Retinex models and robust with noises of different modalities. Various numerical results were implemented to demonstrate the advantages of our model over the existing method for Retinex applications. In the future, we may consider to extend the proposed model to the multiplicative noises, and discuss the convexity and the global convergence of the proposed model.