# A Convex Optimization Model and Algorithm for Retinex

Authors: Qing-Nan Zhao; Ting-Zhu Huang; Xi-Le Zhao; Tian-Hui Ma; Ming-Hui Cheng
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Retinex is a theory on simulating and explaining how human visual system perceives colors under different illumination conditions. The main contribution of this paper is to put forward a new convex optimization model for Retinex. Different from existing methods, the main idea is to rewrite a multiplicative form such that the illumination variable and the reflection variable are decoupled in spatial domain. The resulting objective function involves three terms including the Tikhonov regularization of the illumination component, the total variation regularization of the reciprocal of the reflection component, and the data-fitting term among the input image, the illumination component, and the reciprocal of the reflection component. We develop an alternating direction method of multipliers (ADMM) to solve the convex optimization model. Numerical experiments demonstrate the advantages of the proposed model which can decompose an image into the illumination and the reflection components.


---

## Body

## 1. Introduction

The idea of Retinex was introduced and pioneered by Land and McCann [1] to explain how a combination of processes occurs both in the “retina” and in the “cortex.” Retinex theory tells us that human visual system can ensure that the perceived colors of objects remain relatively constant under varying illumination conditions. That is to say our visual system is robust when it comes to color perception. Retinex theory can deal with the compensation for illumination effect. Therefore, we would like to reduce the influence of nonuniform illumination to enhance an image.Usually, we consider an input image as a two-dimensional functionS which can be decomposed into the illumination function L and the reflection function R. Generally speaking, the input image S is assumed to have the following relation with these two functions:(1)S=L∘R,where ∘ represents the element-wise multiplication.Removing the illumination effect means to decompose the input image into the illumination component and the reflection component. This problem is known to be mathematically ill-posed [2], and many methods have been proposed to solve it in the literature [3–5]. Retinex methods can be classified into random walk methods, recursive methods, center/surround methods, PDE-based methods, and variational methods. Firstly, the original method of Land and McCann was proposed relying on a random walk. The random walk is a discrete time random-process in which the next pixel position is chosen randomly from the neighbors of the current pixel position; see, for example, [6–9]. It needs to regulate many parameters and has high computational complexity. In [10–12], the researchers perform recursive matrix operations to develop recursive methods. The computational efficiency of the recursive methods is improved significantly than that of the random walk methods. However, it is difficult to know how many iterations should be executed in the process. Then Land and McCann put forward the center/surround methods [1]. Later Jobson et al. proposed the single scale Retinex (SSR) and the multiscale Retinex (MSR) [13]. It is easy to implement the SSR and the MSR; both of them need many parameters. In Poisson equation-based methods [14–16], they often convert the original formulation into the logarithmic formulation. Their methods rely on the Mondrian world model which boils down to the assumption on the reflection as a piecewise constant image.Recently, many variational methods were proposed in [17, 18]. The fundamental assumptions are that the illumination component is spatially smooth and the reflection component is piecewise constant. Based on the above assumptions Kimmel et al. presented a variational Retinex formulation [4]. In their model the piecewise constant assumption of the reflection component is not considered. Total Variation (TV) had been widely used in image processing [19–23]. Ma and Osher [24] applied TV and nonlocal TV regularization to Retinex theory. The Bregman iteration was employed to solve their models. It is difficult to set up existence results for their models. Ma et al. further proposed a L1-based variational model to recover the reflection component [25]. In [26], Zosso et al. proposed a unifying framework for Retinex theory. In [27], Liang and Zhang proposed a decomposition model for the Retinex problem via high-order TV.Ng and Wang proposed a TV model for image enhancement [28]. They consider both the illumination component and the reflection component in the objective function. They proposed transforming the multiplicative form (1) into log⁡S=log⁡L+log⁡R. Different from the Ma and Osher model [24], they added some constraints and a fidelity term which ensure that the theoretical analysis can be performed and established. An energy functional was proposed for Retinex as follows:(2)Er,l=∫ΩDr+α2∫Ω∇l2dx+β2∫Ωl-r-s2dx+μ2∫Ωl2dx,where l=log⁡L, r=log⁡R, s=log⁡S, α, β, and μ are positive regularization parameters, and the term ∫Ωl2dx is only used for the theoretical proof which has no practical sense. ∫ΩDr is the total variation term of r [29] and it is employed to characterize the reflection function which makes the model more reasonable. The motivation behind is that we have noticed that once the input image is converted to the logarithmic domain, the small differences between image pixels tend to be ignored. For example, in Figure 1 there is a signal with quite different values in two parts; if we convert the signal to the logarithmic domain, the difference between the two parts is significantly reduced. So we rewrite (1) in spatial domain to avoid the loss of these small textures and details, and at the same time we can also ensure the convexity of the model.Figure 1
(a) a signal; (b) the signal in logarithmic domain.
(a)
(b)The main contribution of this paper is that we develop and study a new convex optimization model for Retinex. A new energy functional is built by considering spatial smoothness of the illumination component and piecewise constant of the reflection component in the decomposition framework. The main idea is to rewrite (1) to a new equation S∘Q=L in spatial domain, where Q is the matrix with [Q]i=1/[R]i. According to the formula after the change, we study this formula in spatial domain which not only makes the illumination variable and the reciprocal of reflection variable decoupled but also can ensure the convexity of the proposed Retinex model, and a convex optimization model can be put forward for Retinex.The rest of this paper is organized as follows. In Section2, we develop a new convex model and employ the numerical method to solve it. In Section 3, we report the experimental results to show the effectiveness of the proposed model and the efficiency of the proposed numerical algorithm. Finally, concluding remarks are given in Section 4.


## 2. The Proposed Model

In this section, we propose a new convex optimization model for Retinex and develop an efficient numerical computational method to solve it. For input RGB color images, we use the HSV color space. The approach is to decompose the color images into three layers: Hue (H), Saturation (S), and Value (V). The V layer corresponds to the brightness of the pixel. Therefore it makes sense to deal with the V layer only since we are trying to extract the illumination image. We denote the V layer of an input imageS by s. s represents the intensity of the light at the pixel. It depends on two main factors. One is the amount of light intensity onto the objects in the image scene, the other is the amount of light source reflected based on the nature of the objects in the image scene. Since s is proportional to the energy radiated by the source, it must be nonzero and finite: 0<s<∞. We consider two functions: the illumination function l and the reflection function r. The V layer of an image s is formed from the two functions:(3)s=l∘r,where we impose the assumptions that 0<l<∞ (illumination effect) and 0<r<1 (reflectivity). The last inequality tells us the reflection image is bounded by 0 (total absorption) and 1 (total reflection). Furthermore, this inequality implies that l>s>0. Based on previous motivation, we rewrite (3) to a new equation (4) in the spatial domain to avoid the loss of these small textures and details(4)s∘q=l,where q is the matrix with [q]i=1/[r]i. Most of existing variational based methods are using the logarithms, which will lose the contrast in both the illumination component l and the reflection component r. Though in the proposed model 1/r will also lose contrast as similar to using the logarithms. However, l and s are both in [0,∞] and keep a good contrast in the proposed model. In addition, the rationality of the proposed model is confirmed by extensive experimental results.Inspired by the convex optimization model for multiplicative noise and blur removal proposed by Zhao et al. [30], a convex optimization model can be put forward for Retinex based on (4).Like all other Retinex algorithms, our Retinex algorithm is based on some basic assumptions as follows:(i) The illumination componentl is spatially smooth. Therefore the Tikhonov regularization term is showed by Dl22, D≜(D(1);D(2)) is the first-order finite difference operator, and the norm Dl2 is defined as Dl2≔∑i,j(D(1)l)i,j2+(D(2)l)i,j2.(ii) Since the reflection componentr is piecewise constant, the reciprocal of the reflection component q is also piecewise constant. Therefore we use the TV to express the regularization term Dq2.(iii) Based on the reflectivity, here we change the constraints into1≤q<∞ and s≤l<∞.A new energy functional is built by considering spatial smoothness of the illumination component and piecewise constant of the reciprocal of the reflection component. In this paper, we propose the following convex optimization energy functional for Retinex to explain how human vision system perceives color. From now on, we will restrict our attention to the discrete setting(5)Eq,l=Dq2+α12Dl22+α22s∘q-l22+δΩ1q+δΩ2l,where α1 and α2 are two positive regularization parameters to control the balance among these terms in the object function, Ω1:{x∣x∈Rn×n,xi,j≥1,∀i,j}; Ω2:{x∣x∈Rn×n,xi,j≥si,j,∀i,j}. Given a nonempty matrix set Ω, the indicator function δΩ(X) is defined by (6)δΩX=0,ifX∈Ω,∞,otherwise. The data fidelity term measured in the l2 norm forces the proximity between s∘q and l. It is reasonable to choose different measures of the fidelity term under different scenarios. For example, the l1 norm can effectively suppress the effect of outliers that may contaminate a given image. This will be considered in our further work.Compared with (2), the new energy functional (5) does not have the last term of (2). The new Retinex algorithm can make the illumination variable and the reciprocal of the reflection variable decoupled, and at the same time it also ensures the convexity of the model. Model (5) is jointly convex for (q,l). Nonconvex models are sensitive to initial values and it is difficult to design and analyze the algorithm. By comparison, our convex model shows the robustness with different initial values and a large number of convex optimization algorithms are available. In addition, all local minimizers are global minimizers and they constitute a closed convex set. For the uniqueness of model (5), although we can not guarantee the uniqueness of the solution, all solutions are corresponding to a global minimum functional value. The associated optimization problem is given by(7)minq,lEq,l.Then we show the existence result for the solution of the proposed model. We first give the following lemma.Lemma 1.
Let{ai} be a sequence of Rn having a bounded gradient; that is, the sequences Dai2 and Dai22 are bounded. If there exists an entry aj that {|aji|} is bounded, then we have that ai2 is bounded.Proof.
For any entryak there exists a set of entries {aj1,aj2,…,ajh} forming a connected road from aj to ak; that is, the entry pairs (aj,aj1), (ajt,ajt+1) for 1≤t≤h-1, and (ajh,ak) are adjacent. Since {ai} has a bounded gradient, then the sequence {Ei}, where(8)Ei=aji-aj1i2+∑t=1h-1ajti-ajt+1i2+ajhi-aki2, is also bounded. Note that there exists a positive constant C satisfying(9)Ei≥Caji-aki2, and then we have that {(aji-aki)2} is bounded, and then {|aki|} is also bounded, which completes the proof.Theorem 2.
Assuming thats∉ Null(D); that is, the null space of D. Model (5) has at least one minimizer.Proof.
It is clear thatE(q,l) is proper and continuous. According to Weierstrass' theorem [31], it remains only to show the coercivity of E(q,l); that is, for every sequence {(qi,li)} such that(10)qi2+li2⟶∞, we have(11)limi→∞⁡Eqi,li=∞. We prove it by contradiction. Suppose that there exists a subsequence of {(qi,li)} (also denoted as {(qi,li)}) that {E(qi,li)} is bounded. According to our assumption that s∉Null(D), there exists two adjacent pixels sj and sj+1 satisfying sj≠sj+1. We consider the following function:(12)E1q,l=2sjqj-lj2+2sj+1qj+1-lj+12. Since {E(qi,li)} is bounded, it is easy to see that the sequences {Dqi2}, {Dli22}, and {E1(qi,li)} are bounded. We have that(13)E1q,l≥sjqj-lj-sj+1qj+1+lj+12=sj-sj+1qj+sj+1qj-qj+1-lj-lj+12. Since {Dqi2} and {Dli22} are bounded, we have that {|qji-qj+1i|} and {|lji-lj+1i|} are bounded. From the assumption sj≠sj+1 we then deduce that {|qji|} is bounded. Moreover, the boundedness of {E1(qi,li)} implies that {|lji|} is also bounded. From Lemma 1 we have {qi2} and {li2} are bounded, which is a contradiction.


### 2.1. The Numerical Method

Many numerical methods are applicable to solve the proposed convex optimization problem including the Bregman method [32], proximal splitting methods [33], primal-dual methods [34–36], and Douglas-Rachford methods [37]. Owing to the convexity of the proposed model, here we just apply and implement the alternating direction method of multipliers (ADMM) by introducing variables w, u, p, and v as follows:(14)minl,w,q,u,p,vw2+α12Du22+α22p-l22+δΩ1v+δΩ2lsubject to w=Dq, l=u, s∘v=p, and v=q.We rewrite (7) to (14), which is the constrained optimization problem. Because all the variables are separated into two groups, (l,w,v) and (u,p,q), we can give the matrix form of the linear constraints: (15)I0000s0I000Ilwv+-I000-I000-D00-Iupq=0, and this form completely accords with the ADMM framework. Here two linear operators both have full column rank, which ensure the convergence of ADMM iteration. We give the resulting augmented Lagrangian function as follows:(16)Ll,v,w,u,p,q,λ=w2+α12Du22+α22p-l22+δΩ1v+δΩ2l+λ1λ2λ3λ4,I0000s0I000Ilwv+-I000-I000-D00-Iupq+β2I0000s0I000Ilwv+-I000-I000-D00-Iupq22,where λ1, λ2, λ3, and λ4 are Lagrange multipliers and β is the penalty parameter [38]. The joint minimization problem can be decomposed into two easier and smaller subproblems such that two groups of variables can be minimized in an alternating order, followed by the update of Lagrange multipliers.InStep1, we update the variables l, w, and v. The optimal solution of the variables can be computed separately since the variables are decoupled. For convenience we omit the superscripts.Thel-subproblem is given by(17)minl⁡α22p-l22+β2l-u+λ1β22+δΩ2l. We can get the solution as follows:(18)lik+1=max⁡α2pk+βuk-λ1k/βiα2+β,si,i=1,2,…,m, where [l]ik+1∈R2;Thew-subproblem is given by(19)minww2+β2w-Dq+λ3β22.We can get the solution as follows:(20)wik+1=max⁡Dqk-λ3kβi2-1β,0Dqk-λ3k/βiDqk-λ3k/βi2,i=1,2,…,m,where wik+1∈R2 and 0·(0/0)=0 is assumed;Thev-subproblem is given by(21)minvβ2s∘v-p+λ2β22+β2v-q+λ4β22+δΩ1v.We can get the solution as follows:(22)vik+1=max⁡s∘pk-λ2k/β+qk-λ4k/βiI+sT∘si,i,1,i=1,2,…,m,where sT denotes the transpose of s and the superscript k denotes the kth iteration of the ADMM.InStep  2, we update the variables u, p, and q. On the same condition, as the variables are decoupled, we can calculate their optimal solutions separately.Theu-subproblem is given by(23)minuα12Du22+β2l-u+λ1β22.We can obtain the corresponding solution with respect to u by solving the normal equation as follows:(24)α1DTD+βu=βl+λ1k,where DTD can be diagonalized by 2D discrete Fourier transforms with periodic boundary conditions. We solve (24) in three steps. First, we compute the right-hand side vector and apply a forward Fast Fourier Transformation (FFT). Second, we get Fuk+1 by dividing the eigenvalues of α1DTD+β, where F denotes the two-dimensional discrete Fourier transform. Third, in order to get the new iterate uk+1, we use the two-dimensional inverse discrete Fourier transform to Fuk+1. At each iteration, there are two FFTs (including one inverse FFT). In this work, we only consider periodic boundary condition by taking advantage of fast Fourier transform (FFT). Figure 2 displays the estimated illumination components with periodic and reflective BCs. We observed that the results are almost the same with periodic and reflective BCs, despite some slight artifacts near boundary with periodic BCs. We believe that more reasonable boundary conditions (BCs) [39, 40] can improve the results.Figure 2
From (a) to (c), the input image, the estimated illumination component with periodic BCs, and the estimated illumination component with reflective BCs.We can get the solution as follows:(25)uk+1=F-1Fβl+λ1kFα1DTD+β.Thep-subproblem is given by(26)minpα22p-l22+β2s∘v-p+λ2β22.We can get the solution as follows:(27)pik+1=α2lk+βs∘vk+λ2kiα2+β,i=1,2,…,m,where [p]ik+1∈R2 is assumed;Theq-subproblem is given by(28)minqβ2w-Dq+λ3β22+β2v-q+λ4β22.As we describe above, we can also obtain the corresponding solution with respect to q by solving the normal equation as follows:(29)DTD+Iq=DTwk+λ3kβ+vk+λ4kβ;here we also use two FFTs to solve this normal equation.We can get the solution as follows:(30)qk+1=F-1FDTwk+λ3k/β+vk+λ4k/βFDTD+I.Finally, we summarize the implementation of ADMM in Algorithm1. The convergence of ADMM is theoretically guaranteed [41, 42]. We use the relative change of the successive iterates as the stopping criterion. The relative change can reflect the convergence behavior and is easy to calculate. In Figure 3 we show the numerical convergence curve to verify the convergence of the proposed ADMM scheme. In the next section, we will show the results of numerical experiments to demonstrate the advantages of the proposed model.Algorithm 1:ADMM procedure for solving (16).
Input. k>0, α1>0, α2>0 and l≥s.
Output. l.
whilestopping criterion lk+1-lk/lk+1≤ϵl is not satisfied do
Step 1. Computing lk+1, wk+1, and vk+1 by solving
m
i
n
l
,
w
,
v
L
(
l
,
v
,
w
,
u
k
,
p
k
,
q
k
,
λ
k
)
.
Step 2. Computing uk+1, pk+1, and qk+1 by solving
m
i
n
u
,
p
,
q
L
(
l
k
+
1
,
v
k
+
1
,
w
k
+
1
,
u
,
p
,
q
,
λ
k
)
.
Step 3. Update λ1k+1, λ2k+1, λ3k+1 and λ4k+1 by the following:
λ
1
k
+
1
⟵
λ
1
k
+
β
(
l
k
+
1
-
u
k
+
1
)
,
λ
2
k
+
1
⟵
λ
2
k
+
β
(
s
∘
v
k
+
1
-
p
k
+
1
)
,
λ
3
k
+
1
⟵
λ
3
k
+
β
(
w
k
+
1
-
D
q
k
+
1
)
,
λ
4
k
+
1
⟵
λ
4
k
+
β
(
v
k
+
1
-
q
k
+
1
)
.
end whileFigure 3
Relative error versus iterations.


## 2.1. The Numerical Method

Many numerical methods are applicable to solve the proposed convex optimization problem including the Bregman method [32], proximal splitting methods [33], primal-dual methods [34–36], and Douglas-Rachford methods [37]. Owing to the convexity of the proposed model, here we just apply and implement the alternating direction method of multipliers (ADMM) by introducing variables w, u, p, and v as follows:(14)minl,w,q,u,p,vw2+α12Du22+α22p-l22+δΩ1v+δΩ2lsubject to w=Dq, l=u, s∘v=p, and v=q.We rewrite (7) to (14), which is the constrained optimization problem. Because all the variables are separated into two groups, (l,w,v) and (u,p,q), we can give the matrix form of the linear constraints: (15)I0000s0I000Ilwv+-I000-I000-D00-Iupq=0, and this form completely accords with the ADMM framework. Here two linear operators both have full column rank, which ensure the convergence of ADMM iteration. We give the resulting augmented Lagrangian function as follows:(16)Ll,v,w,u,p,q,λ=w2+α12Du22+α22p-l22+δΩ1v+δΩ2l+λ1λ2λ3λ4,I0000s0I000Ilwv+-I000-I000-D00-Iupq+β2I0000s0I000Ilwv+-I000-I000-D00-Iupq22,where λ1, λ2, λ3, and λ4 are Lagrange multipliers and β is the penalty parameter [38]. The joint minimization problem can be decomposed into two easier and smaller subproblems such that two groups of variables can be minimized in an alternating order, followed by the update of Lagrange multipliers.InStep1, we update the variables l, w, and v. The optimal solution of the variables can be computed separately since the variables are decoupled. For convenience we omit the superscripts.Thel-subproblem is given by(17)minl⁡α22p-l22+β2l-u+λ1β22+δΩ2l. We can get the solution as follows:(18)lik+1=max⁡α2pk+βuk-λ1k/βiα2+β,si,i=1,2,…,m, where [l]ik+1∈R2;Thew-subproblem is given by(19)minww2+β2w-Dq+λ3β22.We can get the solution as follows:(20)wik+1=max⁡Dqk-λ3kβi2-1β,0Dqk-λ3k/βiDqk-λ3k/βi2,i=1,2,…,m,where wik+1∈R2 and 0·(0/0)=0 is assumed;Thev-subproblem is given by(21)minvβ2s∘v-p+λ2β22+β2v-q+λ4β22+δΩ1v.We can get the solution as follows:(22)vik+1=max⁡s∘pk-λ2k/β+qk-λ4k/βiI+sT∘si,i,1,i=1,2,…,m,where sT denotes the transpose of s and the superscript k denotes the kth iteration of the ADMM.InStep  2, we update the variables u, p, and q. On the same condition, as the variables are decoupled, we can calculate their optimal solutions separately.Theu-subproblem is given by(23)minuα12Du22+β2l-u+λ1β22.We can obtain the corresponding solution with respect to u by solving the normal equation as follows:(24)α1DTD+βu=βl+λ1k,where DTD can be diagonalized by 2D discrete Fourier transforms with periodic boundary conditions. We solve (24) in three steps. First, we compute the right-hand side vector and apply a forward Fast Fourier Transformation (FFT). Second, we get Fuk+1 by dividing the eigenvalues of α1DTD+β, where F denotes the two-dimensional discrete Fourier transform. Third, in order to get the new iterate uk+1, we use the two-dimensional inverse discrete Fourier transform to Fuk+1. At each iteration, there are two FFTs (including one inverse FFT). In this work, we only consider periodic boundary condition by taking advantage of fast Fourier transform (FFT). Figure 2 displays the estimated illumination components with periodic and reflective BCs. We observed that the results are almost the same with periodic and reflective BCs, despite some slight artifacts near boundary with periodic BCs. We believe that more reasonable boundary conditions (BCs) [39, 40] can improve the results.Figure 2
From (a) to (c), the input image, the estimated illumination component with periodic BCs, and the estimated illumination component with reflective BCs.We can get the solution as follows:(25)uk+1=F-1Fβl+λ1kFα1DTD+β.Thep-subproblem is given by(26)minpα22p-l22+β2s∘v-p+λ2β22.We can get the solution as follows:(27)pik+1=α2lk+βs∘vk+λ2kiα2+β,i=1,2,…,m,where [p]ik+1∈R2 is assumed;Theq-subproblem is given by(28)minqβ2w-Dq+λ3β22+β2v-q+λ4β22.As we describe above, we can also obtain the corresponding solution with respect to q by solving the normal equation as follows:(29)DTD+Iq=DTwk+λ3kβ+vk+λ4kβ;here we also use two FFTs to solve this normal equation.We can get the solution as follows:(30)qk+1=F-1FDTwk+λ3k/β+vk+λ4k/βFDTD+I.Finally, we summarize the implementation of ADMM in Algorithm1. The convergence of ADMM is theoretically guaranteed [41, 42]. We use the relative change of the successive iterates as the stopping criterion. The relative change can reflect the convergence behavior and is easy to calculate. In Figure 3 we show the numerical convergence curve to verify the convergence of the proposed ADMM scheme. In the next section, we will show the results of numerical experiments to demonstrate the advantages of the proposed model.Algorithm 1:ADMM procedure for solving (16).
Input. k>0, α1>0, α2>0 and l≥s.
Output. l.
whilestopping criterion lk+1-lk/lk+1≤ϵl is not satisfied do
Step 1. Computing lk+1, wk+1, and vk+1 by solving
m
i
n
l
,
w
,
v
L
(
l
,
v
,
w
,
u
k
,
p
k
,
q
k
,
λ
k
)
.
Step 2. Computing uk+1, pk+1, and qk+1 by solving
m
i
n
u
,
p
,
q
L
(
l
k
+
1
,
v
k
+
1
,
w
k
+
1
,
u
,
p
,
q
,
λ
k
)
.
Step 3. Update λ1k+1, λ2k+1, λ3k+1 and λ4k+1 by the following:
λ
1
k
+
1
⟵
λ
1
k
+
β
(
l
k
+
1
-
u
k
+
1
)
,
λ
2
k
+
1
⟵
λ
2
k
+
β
(
s
∘
v
k
+
1
-
p
k
+
1
)
,
λ
3
k
+
1
⟵
λ
3
k
+
β
(
w
k
+
1
-
D
q
k
+
1
)
,
λ
4
k
+
1
⟵
λ
4
k
+
β
(
v
k
+
1
-
q
k
+
1
)
.
end whileFigure 3
Relative error versus iterations.


## 3. Numerical Experiments

In this section, we present the experiments of Algorithm1 to demonstrate the effectiveness of the proposed Retinex model. For simplicity, we apply our model for color images in HSV color space. We only consider the V channel of HSV color space and then combine the resulting V channel and the original H and S channels back to the resulting color image.As pointed out by [4], we usually get an over-enhanced image by the Retinex process. We can explain the phenomenon by two facts. The first one is that human visual system usually prefers some illumination in the image. The other one is that if we remove all the illumination, some noise may expose in the image from the darker regions of the input image. Therefore, we consider adding a Gamma correction operation on the reconstructed illumination back to the reconstructed reflectance image. From Algorithm 1 we obtain the illumination function l and the initial image s. As discussed in previous sections, the reflectance image r=s/l. The Gamma correction of l with an a parameter γ is performed by l′=w(l/w)1/γ, where w refers to the white value which is equal to 255 in 8 bit images. Therefore, the resulting result s′ is given by(31)s′=l′∘r=l′ls=wl/w1/γls=sl/w1-1/γ,(i)γ=1, the whole illumination is added back, and therefore s′=s.(ii)γ=∞, no illumination is returned, and we get s′=r·w, which is equal to the same reflectance image r, as obtained by the original Retinex. This case is that we add a maximum constant illumination w to the reflectance image r.All experiments are performed in MATLAB R2013a on a PC with a 32-bit operating system and the following configuration: Iter(R) Core(TM) i3-2130 CPU 3.40 GHz and 4 GB RAM. In Figures4(a) and 4(d), we show two input underexposed images for enhancement. In Figures 4(b) and 4(e), the pictures use Gamma correction; we can see the resulting image can not be enhanced. However, in Figures 4(c) and 4(f), we see the pictures using both the proposed Retinex model and Gamma correction can be enhanced to be more natural.Figure 4
(a) and (d) are the input images; (b) and (e) are the resulting images only using Gamma correction; (c) and (f) are the resulting enhanced images using both the proposed Retinex model and Gamma correction.In terms of the human visual system of cognitive function for the image, we give several subjective visual evaluation standards in this paper. At first, the enhanced image should have the appropriate illumination. Then it must have the appropriate color information. Further the appropriate color information should have appropriate spatial distribution which is called the contrast. Finally under the circumstances of undistorted images, we want to have a noise level as lower as possible. In conclusion, we will combine the four visual quality evaluation standard and then make the comprehensive evaluation of the quality for an image.


### 3.1. Experimental Results on Natural Images

In the previous test, the effect of the Gamma correction has been proved. In this subsection, we compare the results by the proposed Retinex model and those by the TV model in [28]. We show five input underexposed images for enhancement. In the tests, we choose the best set of parameters (α1,α2,β,γ)=(1,0.1,10-5,2.2) for the TV model in [28]. For Figure 5, we use (α1,α2,β,γ)=(30,1,200,2.2) for the proposed Retinex model. We set ϵl=1×10-3 for the stopping criteria. The input image is too dark, but the improvement by the proposed Retinex model is clearly visible. The green color of Rubik's cube puzzle is more uniform and natural. It is consistent with the existence of the normal law of things. We repeat the experiments for 20 times and average the results. We find the speeds of the proposed Retinex model are faster than the TV model. Then we show the average results of 20 tests for the iteration number and CPU time (s) in Table 1.Table 1
The iteration number and CPU time (s) of the TV model and the proposed Retinex model.
The TV model
The proposed Retinex model
Iteration number
39
97
Time (s)
4.92
3.13Figure 5
The first row from left to right: the input image, enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), and the proposed Retinex model with the parameters (α1,α2,β,γ)=(30,1,200,2.2); the second row: zoom into the enhanced images.In Figure6, the influence of illumination is removed by the TV model and the proposed Retinex model as can be seen in the first row. However, zooming into the images we can see that more shadows are removed by the proposed Retinex model with the parameters (α1,α2,β,γ)=(10,1,200,2.2). The letters hidden in the shadows are clearly visible. The zoom results show that the color channels are recovered very well in the whole image. Figure 7 shows another image enhancement example; the proposed Retinex model with parameters (α1,α2,β,γ)=(20,1,200,2.2) provides a smoother image without destroying important image details. Zooming into the images, the lamps in the wall are brighter without the shadows on the top of the image.Figure 6
The first row from left to right: the input image, enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), and the proposed Retinex model with the parameters (α1,α2,β,γ)=(10,1,200,2.2); the second row: zoom into the enhanced images.Figure 7
The first row from left to right: the input image, enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), and the proposed Retinex model with the parameters (α1,α2,β,γ)=(20,1,200,2.2); the second row: zoom into the enhanced images.Figure8 shows an example where the illumination is corrected. The proposed Retinex model with the parameters (α1,α2,β,γ)=(10,2,200,2.2) is used to separate the illumination and the reflection. The output enhanced images are sharper and clearer than those by the TV method. Zooming into the images, we can see the colors of the maple leaves are better than those using the TV model. Finally, in Figure 9 our experimental results are performed by the proposed Retinex model very well. The visual effects of the enhanced images by the proposed Retinex model with the parameters (α1,α2,β,γ)=(20,3,200,2.2) are better than those by the TV method. The zoom results also make this point clear; the windows are clearer and brighter. We also show the comparisons of CPU time (s) in Table 2.Table 2
The CPU time (s) of the TV model and the proposed Retinex model.
The TV model
The proposed Retinex model
Image Book
11.37
2.63
Image Room
17.16
10.28
Image God
3.05
1.65
Image House
41.51
30.67Figure 8
The first row from left to right: the input image, enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), and the proposed Retinex model with the parameters (α1,α2,β,γ)=(10,2,200,2.2); the second row: zoom into the enhanced images.Figure 9
The first row from left to right: the input image, enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), and the proposed Retinex model with the parameters (α1,α2,β,γ)=(20,3,200,2.2); the second row: zoom into the enhanced images.


### 3.2. The Effect of Parameters

In the previous test, we compare the results by the proposed Retinex model and those by the TV model. The enhanced images have good visual effects using the proposed Retinex model. In Figure10, we test the proposed Retinex model using different values of γ: 2.2, 10, 100 to get the enhanced images. We use the best parameters of TV model with (α1,α2,β,γ)=(1,0.1,10-5,2.2). In the proposed Retinex model, we use the parameters (α1,α2,β,γ)=(10,10,10,2.2). We find that the enhanced images using the TV model are not enough sharp and pleasant. Particularly when γ is large, some details such as the side of the street are lost. By comparison, the enhanced images using the proposed Retinex model are visually appealing, and the proposed Retinex method is more robust for different values of γ.Figure 10
The first row: the input image and the enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), respectively; the second row: the input image and the enhanced images by the proposed Retinex model with (α1,α2,β,γ)=(10,10,10,2.2), respectively.In the next test, we discuss how to choose the parameters for the proposed Retinex model. Then we study the effect of parametersα1 and α2 of the proposed Retinex model on the results of the experiments. At first, we set α2=10 and α1=0.1,1,10,100, respectively. In the first row of Figure 11, when the parameter α1 increased, the reflection component becomes sharper, and the second row shows that the illumination component becomes smoother. Then we set α1=10 and α2=1,10,100,1000, respectively. In the first row of Figure 12, when the parameter α2 increased, the reflection component becomes flatter, and the second row shows that the illumination component becomes more rough.Figure 11
The first row: the reflection component by the proposed Retinex model withα2=10 and α1=0.1,1,10,100, respectively; the second row: the illumination component by the proposed Retinex model with α2=10 and α1=0.1,1,10,100, respectively.Figure 12
The first row: the reflection component by the proposed Retinex model withα1=10 and α2=1,10,100,1000, respectively; the second row: the illumination component by the proposed Retinex model with α1=10 and α2=1,10,100,1000, respectively.


### 3.3. Experimental Results on Medical Images

In the same way, we test some medical images as well in Figure13 and compare the proposed Retinex model with the TV model in [28]. In particular, we can see more details in the brain from the results by the proposed Retinex model than those by the TV model. Both gray and white matter are clear. Skeletal outline is clearly visible. Therefore, the proposed Retinex model can generate better enhanced images with high contrast and preserve more features in the enhancement. Numerical experiments demonstrate the effectiveness of the proposed model and the efficiency of the proposed numerical algorithm. We show the comparisons of CPU time (s) in Table 3 as well.Table 3
The CPU time (s) of the TV model and the proposed Retinex model.
The TV model
The proposed Retinex model
Image Brain
4.19
2.09
Image Bones
24.75
5.67Figure 13
(a) and (d) are the input images; (b) and (e) are the resulting images by the TV model; (c) and (f) are the resulting images by the proposed Retinex model.


## 3.1. Experimental Results on Natural Images

In the previous test, the effect of the Gamma correction has been proved. In this subsection, we compare the results by the proposed Retinex model and those by the TV model in [28]. We show five input underexposed images for enhancement. In the tests, we choose the best set of parameters (α1,α2,β,γ)=(1,0.1,10-5,2.2) for the TV model in [28]. For Figure 5, we use (α1,α2,β,γ)=(30,1,200,2.2) for the proposed Retinex model. We set ϵl=1×10-3 for the stopping criteria. The input image is too dark, but the improvement by the proposed Retinex model is clearly visible. The green color of Rubik's cube puzzle is more uniform and natural. It is consistent with the existence of the normal law of things. We repeat the experiments for 20 times and average the results. We find the speeds of the proposed Retinex model are faster than the TV model. Then we show the average results of 20 tests for the iteration number and CPU time (s) in Table 1.Table 1
The iteration number and CPU time (s) of the TV model and the proposed Retinex model.
The TV model
The proposed Retinex model
Iteration number
39
97
Time (s)
4.92
3.13Figure 5
The first row from left to right: the input image, enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), and the proposed Retinex model with the parameters (α1,α2,β,γ)=(30,1,200,2.2); the second row: zoom into the enhanced images.In Figure6, the influence of illumination is removed by the TV model and the proposed Retinex model as can be seen in the first row. However, zooming into the images we can see that more shadows are removed by the proposed Retinex model with the parameters (α1,α2,β,γ)=(10,1,200,2.2). The letters hidden in the shadows are clearly visible. The zoom results show that the color channels are recovered very well in the whole image. Figure 7 shows another image enhancement example; the proposed Retinex model with parameters (α1,α2,β,γ)=(20,1,200,2.2) provides a smoother image without destroying important image details. Zooming into the images, the lamps in the wall are brighter without the shadows on the top of the image.Figure 6
The first row from left to right: the input image, enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), and the proposed Retinex model with the parameters (α1,α2,β,γ)=(10,1,200,2.2); the second row: zoom into the enhanced images.Figure 7
The first row from left to right: the input image, enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), and the proposed Retinex model with the parameters (α1,α2,β,γ)=(20,1,200,2.2); the second row: zoom into the enhanced images.Figure8 shows an example where the illumination is corrected. The proposed Retinex model with the parameters (α1,α2,β,γ)=(10,2,200,2.2) is used to separate the illumination and the reflection. The output enhanced images are sharper and clearer than those by the TV method. Zooming into the images, we can see the colors of the maple leaves are better than those using the TV model. Finally, in Figure 9 our experimental results are performed by the proposed Retinex model very well. The visual effects of the enhanced images by the proposed Retinex model with the parameters (α1,α2,β,γ)=(20,3,200,2.2) are better than those by the TV method. The zoom results also make this point clear; the windows are clearer and brighter. We also show the comparisons of CPU time (s) in Table 2.Table 2
The CPU time (s) of the TV model and the proposed Retinex model.
The TV model
The proposed Retinex model
Image Book
11.37
2.63
Image Room
17.16
10.28
Image God
3.05
1.65
Image House
41.51
30.67Figure 8
The first row from left to right: the input image, enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), and the proposed Retinex model with the parameters (α1,α2,β,γ)=(10,2,200,2.2); the second row: zoom into the enhanced images.Figure 9
The first row from left to right: the input image, enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), and the proposed Retinex model with the parameters (α1,α2,β,γ)=(20,3,200,2.2); the second row: zoom into the enhanced images.


## 3.2. The Effect of Parameters

In the previous test, we compare the results by the proposed Retinex model and those by the TV model. The enhanced images have good visual effects using the proposed Retinex model. In Figure10, we test the proposed Retinex model using different values of γ: 2.2, 10, 100 to get the enhanced images. We use the best parameters of TV model with (α1,α2,β,γ)=(1,0.1,10-5,2.2). In the proposed Retinex model, we use the parameters (α1,α2,β,γ)=(10,10,10,2.2). We find that the enhanced images using the TV model are not enough sharp and pleasant. Particularly when γ is large, some details such as the side of the street are lost. By comparison, the enhanced images using the proposed Retinex model are visually appealing, and the proposed Retinex method is more robust for different values of γ.Figure 10
The first row: the input image and the enhanced images by the TV model with the best set of parameters(α1,α2,β,γ)=(1,0.1,10-5,2.2), respectively; the second row: the input image and the enhanced images by the proposed Retinex model with (α1,α2,β,γ)=(10,10,10,2.2), respectively.In the next test, we discuss how to choose the parameters for the proposed Retinex model. Then we study the effect of parametersα1 and α2 of the proposed Retinex model on the results of the experiments. At first, we set α2=10 and α1=0.1,1,10,100, respectively. In the first row of Figure 11, when the parameter α1 increased, the reflection component becomes sharper, and the second row shows that the illumination component becomes smoother. Then we set α1=10 and α2=1,10,100,1000, respectively. In the first row of Figure 12, when the parameter α2 increased, the reflection component becomes flatter, and the second row shows that the illumination component becomes more rough.Figure 11
The first row: the reflection component by the proposed Retinex model withα2=10 and α1=0.1,1,10,100, respectively; the second row: the illumination component by the proposed Retinex model with α2=10 and α1=0.1,1,10,100, respectively.Figure 12
The first row: the reflection component by the proposed Retinex model withα1=10 and α2=1,10,100,1000, respectively; the second row: the illumination component by the proposed Retinex model with α1=10 and α2=1,10,100,1000, respectively.


## 3.3. Experimental Results on Medical Images

In the same way, we test some medical images as well in Figure13 and compare the proposed Retinex model with the TV model in [28]. In particular, we can see more details in the brain from the results by the proposed Retinex model than those by the TV model. Both gray and white matter are clear. Skeletal outline is clearly visible. Therefore, the proposed Retinex model can generate better enhanced images with high contrast and preserve more features in the enhancement. Numerical experiments demonstrate the effectiveness of the proposed model and the efficiency of the proposed numerical algorithm. We show the comparisons of CPU time (s) in Table 3 as well.Table 3
The CPU time (s) of the TV model and the proposed Retinex model.
The TV model
The proposed Retinex model
Image Brain
4.19
2.09
Image Bones
24.75
5.67Figure 13
(a) and (d) are the input images; (b) and (e) are the resulting images by the TV model; (c) and (f) are the resulting images by the proposed Retinex model.


## 4. Conclusion

In this paper, we presented a new convex optimization model for Retinex. Different from the existing methods, we developed a new framework so as to make the illumination variable and the reflection variable decoupled. We showed the existence of the solution of the proposed Retinex model. In particular, we employed a computation method ADMM to solve the proposed Retinex model. Numerical experimental results were given to illustrate that the proposed model and the computation method are effective and efficient to improve the quality of the enhanced images.



---

*Source: 4012767-2017-07-24.xml*