A Convex Optimization Model and Algorithm for Retinex

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.


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 function  which can be decomposed into the illumination function  and the reflection function .Generally speaking, the input image  is assumed to have the following relation with these two functions: 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][4][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][7][8][9].It needs to regulate many parameters and has high computational complexity.In [10][11][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][15][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][20][21][22][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  1 -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  = log  + log .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: where  = log ,  = log ,  = log , , , and  are positive regularization parameters, and the term ∫ Ω  2  is only used for the theoretical proof which has no practical sense.∫ Ω || is the total variation term of  [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.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  ∘  =  in spatial domain, where  is the matrix with []  = 1/[]  .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 Section 2, 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.

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 image  by . 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  is proportional to the energy radiated by the source, it must be nonzero and finite: 0 <  < ∞.We consider two functions: the illumination function  and the reflection function .The V layer of an image  is formed from the two functions: where we impose the assumptions that 0 <  < ∞ (illumination effect) and 0 <  < 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  >  > 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 where  is the matrix with []  = 1/[]  .Most of existing variational based methods are using the logarithms, which will lose the contrast in both the illumination component  and the reflection component .Though in the proposed model 1/ will also lose contrast as similar to using the logarithms.However,  and  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 component  is spatially smooth.Therefore the Tikhonov regularization term is showed by ‖‖ 2  2 ,  ≜ ( (1) ;  (2) ) is the first-order finite difference operator, and the norm ‖‖ 2 is defined as ‖‖ 2 fl ∑ , √( (1) ) 2  , + ( (2) ) 2 , .(ii) Since the reflection component  is piecewise constant, the reciprocal of the reflection component  is also piecewise constant.Therefore we use the TV to express the regularization term ‖‖ 2 .(iii) Based on the reflectivity, here we change the constraints into 1 ≤  < ∞ and  ≤  < ∞.
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 where  1 and  2 are two positive regularization parameters to control the balance among these terms in the object function, The data fidelity term measured in the ℓ 2 norm forces the proximity between  ∘  and .It is reasonable to choose different measures of the fidelity term under different scenarios.For example, the ℓ 1 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 (, ).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 min ,  (, ) . ( Then we show the existence result for the solution of the proposed model.We first give the following lemma.
is also bounded.Note that there exists a positive constant  satisfying and then we have that {(   −    ) 2 } is bounded, and then {|   |} is also bounded, which completes the proof.Theorem 2. Assuming that  ∉ Null(); that is, the null space of .Model ( 5) has at least one minimizer.
We rewrite (7) to (14), which is the constrained optimization problem.Because all the variables are separated into two groups, (, , V) and (, , ), we can give the matrix form of the linear constraints: (15) 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: 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.
In Step 1, we update the variables , , and V.The optimal solution of the variables can be computed separately since the variables are decoupled.For convenience we omit the superscripts. The We can get the solution as follows: where The -subproblem is given by We can get the solution as follows: [] We can get the solution as follows: where   denotes the transpose of  and the superscript  denotes the th iteration of the ADMM.
In Step 2, we update the variables , , and .On the same condition, as the variables are decoupled, we can calculate their optimal solutions separately.
The -subproblem is given by We can obtain the corresponding solution with respect to  by solving the normal equation as follows: where    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 F +1 by dividing the eigenvalues of  1    + , where F denotes the two-dimensional discrete Fourier transform.Third, in order to get the new iterate  +1 , we use the two-dimensional inverse discrete Fourier transform to F +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.We can get the solution as follows: The -subproblem is given by We can get the solution as follows: [] where As we describe above, we can also obtain the corresponding solution with respect to  by solving the normal equation as follows: here we also use two FFTs to solve this normal equation.We can get the solution as follows: ) . ( Finally, we summarize the implementation of ADMM in Algorithm 1.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.

Numerical Experiments
In this section, we present the experiments of Algorithm 1 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  Input. > 0,  1 > 0,  2 > 0 and  ≥ .Output..
Step 3. Update  +1 1 ,  +1 2 ,  +1 3 and  +1 4 by the following: ). end while Algorithm 1: ADMM procedure for solving (16).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  and the initial image .As discussed in previous sections, the reflectance image  = /.The Gamma correction of  with an a parameter  is performed by   = (/) 1/ , where  refers to the white value which is equal to 255 in 8 bit images.Therefore, the resulting result   is given by (i)  = 1, the whole illumination is added back, and therefore   = .
(ii)  = ∞, no illumination is returned, and we get   = ⋅ , which is equal to the same reflectance image , as obtained by the original Retinex.This case is that we add a maximum constant illumination  to the reflectance image .
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 Figures 4(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.
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.

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   = 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.
In Figure 6, 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 8 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.

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 Figure 10, 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 .
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.

Experimental Results on Medical
Images.In the same way, we test some medical images as well in Figure 13 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.

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.

Figure 1 :
Figure 1: (a) a signal; (b) the signal in logarithmic domain.

Figure 2 :Figure 3 :
Figure 2: From (a) to (c), the input image, the estimated illumination component with periodic BCs, and the estimated illumination component with reflective BCs.

Figure 4 :
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.

Figure 13 :
Figure 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.

Table 1 :
The iteration number and CPU time (s) of the TV model and the proposed Retinex model.

Table 2 :
The CPU time (s) of the TV model and the proposed Retinex model.

Table 3 :
The CPU time (s) of the TV model and the proposed Retinex model.