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.

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 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 [35]. 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, [69]. It needs to regulate many parameters and has high computational complexity. In [1012], 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 [1416], 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 [1923]. 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 -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 . 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 , , , , , and are positive regularization parameters, and the term 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 . 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.

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 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: . 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 (illumination effect) and (reflectivity). The last inequality tells us the reflection image is bounded by 0 (total absorption) and 1 (total reflection). Furthermore, this inequality implies that . 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 detailswhere is the matrix with . 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 will also lose contrast as similar to using the logarithms. However, and are both in 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 , is the first-order finite difference operator, and the norm is defined as .

(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 .

(iii) Based on the reflectivity, here we change the constraints into 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 settingwhere and are two positive regularization parameters to control the balance among these terms in the object function, ; . Given a nonempty matrix set , the indicator function is defined by The data fidelity term measured in the norm forces the proximity between and . It is reasonable to choose different measures of the fidelity term under different scenarios. For example, the 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 byThen we show the existence result for the solution of the proposed model. We first give the following lemma.

Lemma 1. Let be a sequence of having a bounded gradient; that is, the sequences and are bounded. If there exists an entry that is bounded, then we have that is bounded.

Proof. For any entry there exists a set of entries forming a connected road from to ; that is, the entry pairs , for , and are adjacent. Since has a bounded gradient, then the sequence , where is also bounded. Note that there exists a positive constant satisfying and then we have that 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.

Proof. It is clear that is proper and continuous. According to Weierstrass' theorem [31], it remains only to show the coercivity of ; that is, for every sequence such that we have We prove it by contradiction. Suppose that there exists a subsequence of (also denoted as ) that is bounded. According to our assumption that , there exists two adjacent pixels and satisfying . We consider the following function: Since is bounded, it is easy to see that the sequences , , and are bounded. We have that Since and are bounded, we have that and are bounded. From the assumption we then deduce that is bounded. Moreover, the boundedness of implies that is also bounded. From Lemma 1 we have and 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 [3436], 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 , , , and as follows:subject to , , , and .

We rewrite (7) to (14), which is the constrained optimization problem. Because all the variables are separated into two groups, and , we can give the matrix form of the linear constraints: 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 , , , and 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 . The optimal solution of the variables can be computed separately since the variables are decoupled. For convenience we omit the superscripts.

The -subproblem is given by We can get the solution as follows: where ;

The -subproblem is given byWe can get the solution as follows:where and is assumed;

The -subproblem is given byWe 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 byWe 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 by dividing the eigenvalues of , where denotes the two-dimensional discrete Fourier transform. Third, in order to get the new iterate , we use the two-dimensional inverse discrete Fourier transform to . 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 byWe can get the solution as follows:where is assumed;

The -subproblem is given byAs 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.

Input. , , and .
Output. .
while stopping criterion is not satisfied do
Step 1. Computing , , and by solving
Step 2. Computing , , and by solving
Step 3. Update , , and by the following:
end while

3. 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 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 , where refers to the white value which is equal to 255 in 8 bit images. Therefore, the resulting result is given by

(i) , 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.

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 for the TV model in [28]. For Figure 5, we use for the proposed Retinex model. We set 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 . 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 ( 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 ( 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 ( 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.

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 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 . In the proposed Retinex model, we use the parameters . 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 and of the proposed Retinex model on the results of the experiments. At first, we set and , respectively. In the first row of Figure 11, when the parameter increased, the reflection component becomes sharper, and the second row shows that the illumination component becomes smoother. Then we set and , respectively. In the first row of Figure 12, when the parameter increased, the reflection component becomes flatter, and the second row shows that the illumination component becomes more rough.

3.3. 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.

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.

Disclosure

This work is partially included in the first author’s master thesis [43].

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors would like to express their great thankfulness to Dr. Wei Wang, Professor M. K. Ng, and Dr. Huibin Chang for allowing the use of their codes and images. This research is supported by 973 Program (2013CB329404), NSFC (61370147, 61402082), and the Fundamental Research Funds for the Central Universities (ZYGX2016J132).