Realistic Image Rendition Using a Variable Exponent Functional Model for Retinex

The goal of realistic image rendition is to recover the acquired image under imperfect illuminant conditions, where non–uniform illumination may degrade image quality with high contrast and low SNR. In this paper, the assumption regarding illumination is modified and a variable exponent functional model for Retinex is proposed to remove non–uniform illumination and reduce halo artifacts. The theoretical derivation is provided and experimental results are presented to illustrate the effectiveness of the proposed model.


Introduction
Realistic image rendition aims to represent human perception of natural scenes. The meaning of "realistic" is to provide machine vision with ideal images according to the human visual system. A complete visual pathway includes the optic nerve, retina, optic tract, optic chiasm, superior colliculus, lateral geniculate nucleus, optic radiation, and cortex, as shown diagrammatically in Figure 1 [1]. The main features of realistic image rendition include color constancy, image enhancement, high dynamic range compression, etc. The physiological basis for color constancy involves specialized neurons in the primary visual cortex that compute local ratios of cone activity [2], which is the same calculation as Land's Retinex algorithm [3,4] used to achieve color constancy. The existence of these specialized cells, double-opponent cells, has been proven using receptive field mapping. Receptive field [5,6] is the basic unit of visual information processing, and can be separated into two types: On-Center and Off-Center ganglion cells. Figure 2 shows the receptive field in the retina.
Algorithms of realistic image rendition based on visual characteristics generally include Retinex algorithms for color constancy. The word "Retinex" is a combination of "retina" and "cortex". The aim of Retinex theory is to tell whether human eyes can determine reflectance when both the illumination and reflectance are unknown. Land and McCann [6] first proposed the Retinex theory, a path-based algorithm, as a model of color perception of the human visual system (HSV). Many algorithms [7][8][9] are based on this approach, which differ in how the path is selected. However, these methods have high computation complexity and require numerous parameters. McCann [10][11][12] replaced the path calculation by a recursive matrix computation which greatly improved computational efficiency. However, the terminal criterion is not clear and can strongly influence the result. In PDE based models [13], the Retinex principles are often translated into a physical form. These algorithms are developed based on solving a Poisson equation which can yield fast and exact implementation using only two fast Fourier transforms. The main assumption in this algorithm type is that the is that the reflectance performs as the sharp details in the image, while illumination varies smoothly. Based on the assumptions used in PDE formulations, Kimmel et al. [14] proposed a general variational model for the Retinex problem that unified previous methods. Ma and Osher [15,16] proposed a total variation and nonlocal total variation(TV) regularized model using the same assumptions. Ng el at. [17] investigated the TV model with more constraints. Recently, Liang and Zhang [18] established a new higher order total variation L1 decomposition model (HoTVL1) which can correct the piecewise linear shadows. Zosso [19,20] proposed a unifying Retinex model based on non-local differential operators.  is that the reflectance performs as the sharp details in the image, while illumination varies smoothly. Based on the assumptions used in PDE formulations, Kimmel et al. [14] proposed a general variational model for the Retinex problem that unified previous methods. Ma and Osher [15,16] proposed a total variation and nonlocal total variation(TV) regularized model using the same assumptions. Ng el at. [17] investigated the TV model with more constraints. Recently, Liang and Zhang [18] established a new higher order total variation L1 decomposition model (HoTVL1) which can correct the piecewise linear shadows. Zosso [19,20] proposed a unifying Retinex model based on non-local differential operators.   To the best of our knowledge, almost all of the important assumptions about illumination in existing Retinex models require spatial smoothness. However, many images with non-uniform illuminations have non smooth illumination, actually. In this paper, we assume: (a) The reflecting object a Lambertian reflector and reflectance corresponds to sharp details in the image; (b) Illumination is smooth in most regions, but may contain non-smooth part(s).
Based on these assumptions, we propose a new Retinex model using a variable exponent functional. We assume that the illumination function belongs to some Sobolev space with variable exponents. The proposed model solution existence is proved here. Although the proposed model is developed for specific cases, it can also be applied to general degraded images and significantly reduces the halo artifact.
In Section 2, we argue the reasonability of the assumption and present the proposed model. We also present a proof of solution existence for the proposed model, and introduce an efficient iterative solution method. In Section 3, we present several numerical examples to demonstrate the effectiveness of the proposed model. Concluding remarks are presented in Section 4.

New Assumption
To illustrate the proposed assumption, let us consider the images with different illumination conditions and their corresponding surfaces in Figure 3. The corresponding surfaces illustrate the shadow shapes. The illumination of the text image varies smoothly, whereas that of the book image has an apparent non-smooth component. Figure 4 shows a single row from the two images, illustrating the text image curve changes relatively smooth in the shadow area, while that of the book image changes dramatically at the edge of the shadow and relatively smooth in the interior of the shadow. To the best of our knowledge, almost all of the important assumptions about illumination in existing Retinex models require spatial smoothness. However, many images with non-uniform illuminations have non smooth illumination, actually. In this paper, we assume: (a) The reflecting object a Lambertian reflector and reflectance corresponds to sharp details in the image; (b) Illumination is smooth in most regions, but may contain non-smooth part(s).
Based on these assumptions, we propose a new Retinex model using a variable exponent functional. We assume that the illumination function belongs to some Sobolev space with variable exponents. The proposed model solution existence is proved here. Although the proposed model is developed for specific cases, it can also be applied to general degraded images and significantly reduces the halo artifact.
In Section 2, we argue the reasonability of the assumption and present the proposed model. We also present a proof of solution existence for the proposed model, and introduce an efficient iterative solution method. In Section 3, we present several numerical examples to demonstrate the effectiveness of the proposed model. Concluding remarks are presented in Section 4.

New Assumption
To illustrate the proposed assumption, let us consider the images with different illumination conditions and their corresponding surfaces in Figure 3. The corresponding surfaces illustrate the shadow shapes. The illumination of the text image varies smoothly, whereas that of the book image has an apparent non-smooth component. Figure 4 shows a single row from the two images, illustrating the text image curve changes relatively smooth in the shadow area, while that of the book image changes dramatically at the edge of the shadow and relatively smooth in the interior of the shadow.    The above examples support the proposed assumption. Indeed, every severe non-uniform illumination case is likely to have a non-smooth part. Our aim is to extract illumination images and recover realistic images.

Proposed Model
First, we introduce the variable exponent functional and some related models. Blomgren et al. [21] proposed the variable exponent functional for image denoising problems. They tried to minimize: where u is the image function and p is a monotonically decreasing function with lims∞p(s) = 1. Choosing p = 1 produces the widely used Rudin-Osher-Fatemi (ROF) model [22] which preserves edge sharpness, but often causes the "staircase" effect; p = 2 produces isotropic diffusion, which avoids the "staircase" effect but smears edges. Thus, it is natural to combine their benefits with a variable exponent. However, because p relies on u  , it is difficult to establish the lower semi continuity property of the functional. Chenet al. [23] proposed a variable exponent linear growth functional model for image denoising, enhancement and restoration, which is extended by Liet al. [24], using variable exponent functionals in image restoration. For simplicity, we formulate and discuss our model based on grayscale images. For color images, we simply map the color into HSV(hue, saturation, value) color space, process only the V channel, then transform it back to the RGB domain. This method is called HSV Retinex [14,17].
Let I be an image defined in image domain  . The primary goal of Retinex theory is to decompose I into the reflectance image, R, and the illumination image, L, as shown in Figure 5, such that, at each point in the image domain [25]: and following [14,17], we may further assume that: We first convert Equation (2) into the logarithmic domain: The above examples support the proposed assumption. Indeed, every severe non-uniform illumination case is likely to have a non-smooth part. Our aim is to extract illumination images and recover realistic images.

Proposed Model
First, we introduce the variable exponent functional and some related models. Blomgren et al. [21] proposed the variable exponent functional for image denoising problems. They tried to minimize: where u is the image function and p is a monotonically decreasing function with lim sÑ0 ppsq " 2, lim sÑ8 p(s) = 1. Choosing p = 1 produces the widely used Rudin-Osher-Fatemi (ROF) model [22] which preserves edge sharpness, but often causes the "staircase" effect; p = 2 produces isotropic diffusion, which avoids the "staircase" effect but smears edges. Thus, it is natural to combine their benefits with a variable exponent. However, because p relies on ∇u, it is difficult to establish the lower semi continuity property of the functional. Chen et al. [23] proposed a variable exponent linear growth functional model for image denoising, enhancement and restoration, which is extended by Li et al. [24], using variable exponent functionals in image restoration. For simplicity, we formulate and discuss our model based on grayscale images. For color images, we simply map the color into HSV(hue, saturation, value) color space, process only the V channel, then transform it back to the RGB domain. This method is called HSV Retinex [14,17].
Let I be an image defined in image domain Ω. The primary goal of Retinex theory is to decompose I into the reflectance image, R, and the illumination image, L, as shown in Figure 5, such that, at each point in the image domain [25]: and following [14,17], we may further assume that: We first convert Equation (2) into the logarithmic domain: so that: i " l`r Based on our new assumption, the illumination image may contain non-smooth parts. Weuse a total variation like regularizer near non-smooth parts and a Tikhonov like regularizer for smooth parts. We minimize the objective function as follows: where λ is a positive number, and ppxq " 1`1 1`w|∇d| 2 , where d is the ideal illumination image, discussed in Section 2.4.2. The first fidelity term on the right side of model (3) measures the similarity of the gradient between the illumination and the original image, and the second is the regularization term. Clearly, p Ñ 1 near the edges of d where the gradient is large, and so the regularizer is similar to a TV regularizer which can preserve edges; p Ñ 2 in the homogeneous regions where the gradient is small, and here the regularizer is similar to a Tikhonov like regularizer, which is superior to total variation. In other regions, the penalty is adjusted by p(x). Based on our new assumption, the illumination image may contain non-smooth parts. Weuse a total variation like regularizer near non-smooth parts and a Tikhonov like regularizer for smooth parts. We minimize the objective function as follows: where  is a positive number, and The first fidelity term on the right side of model (3) measures the similarity of the gradient between the illumination and the original image, and the second is the regularization term. Clearly, 1 p  near the edges of dwhere the gradient is large, and so the regularizer is similar to a TV regularizer which can preserve edges; 2 p  in the homogeneous regions where the gradient is small, and here the regularizer is similar to a Tikhonov like regularizer, which is superior to total variation. In other regions, the penalty is adjusted by p(x). The classical Retinex algorithm uses a Gaussian filter, equivalent to a Tikhonov regularizer, to obtain the illumination image. However, a Gaussian filter smears edges, which is the main cause of halo artifacts [26]. Using the adaptive TV like regularizer for the high contrast edges in the image, our model not only prevents halo artifacts but also extracts the edges of non-uniform illumination from the image.

Solution Existence
Let us recall some definitions and basic properties of variable exponent Lebesgue and Sobolev spaces, following [24,27]. We define a functional, which is also called modular: and a norm: Then the variable exponent Lebesgue and Sobolev spaces are, respectively: The classical Retinex algorithm uses a Gaussian filter, equivalent to a Tikhonov regularizer, to obtain the illumination image. However, a Gaussian filter smears edges, which is the main cause of halo artifacts [26]. Using the adaptive TV like regularizer for the high contrast edges in the image, our model not only prevents halo artifacts but also extracts the edges of non-uniform illumination from the image.

Solution Existence
Let us recall some definitions and basic properties of variable exponent Lebesgue and Sobolev spaces, following [24,27].

Theorem 1.
Let Ω Ă R 2 be a bounded open set with Lipchitz boundary, i P W 1,ppxq pΩq X L 2 pΩq, then the minimization problem: has a minimizer l P W 1,ppxq pΩq X L 2 pΩq.
Proof of Theorem 1. Let tl k u 8 k"1 , l k P W 1,ppxq pΩq X L 2 pΩq be the minimizing sequence for Eplq. Then: ż Ω |∇l k´∇ i| 2 dx ă M and where M denotes a universal positive constant that may differ from line to line. Hence ş Ω |∇l k | 2 dx ă M.
Therefore, ş Ω |l k | ppxq dx ă M, and together with the inequality ş Ω |∇l k | ppxq dx ă M, we obtain Q ppxq pl k q`Q ppxq p∇l k q ă M. This implies that tl k u 8 k"1 is a uniformly bounded sequence in W 1,ppxq pΩq due to Lemma 1, and tl k u 8 k"1 is also uniformly bounded in L 2 pΩq. Since W 1,ppxq pΩq X L 2 pΩq is a reflexive Banach space, up to a subsequence, there exists l˚P W 1,ppxq pΩq X L 2 pΩq such that ! l k j ) converges weakly to l˚in W 1,ppxq pΩq X L 2 pΩq. From Lemma 4, Eplq is lower semi continuous in W 1,ppxq pΩq X L 2 pΩq. Thus: Epl˚q ď lim kÑ8 Epl k q " inf lPW 1,ppxq pΩqXL 2 pΩq Eplq Therefore, l˚is the minimum point of Eplq.

Implementation
We formulate the basic procedure for solving problem Equation (3) following the split Bregman [28][29][30][31][32][33] technique. We solve the minimization by introducing an auxiliary variable b: By adding one quadratic penalty function term, we convert Equation (4) to an unconstrained splitting formulation: where γ is a positive parameter which controls the weight of the penalty term. Similar to the split Bregman iteration, we propose the scheme: Alternatively, this joint minimization problem can be solved by decomposing into several subproblems.

Subproblem l with Fixed b and t
Given the fixed variable b k and t k , our aim is to find the solution of the problem: which has the optimality condition: pγ`1q∆l " γ∇¨pb k´tk q`∆i (8) where b " pb x , b y q and t " pt x , t y q. Since the discrete system is strictly diagonally dominant with Neumann boundary condition, the most natural choice is the Gauss-Seidel method. The Gauss-Seidel solution to this subproblem can be written componentwise as: Note that this subproblem can also be solved by FFT with periodic boundary condition.

Subproblem b with Fixed l and t
Similarly, we solve: which has the optimality condition: # λppxq|b| ppxq´2 b x`2 γpb x´∇x l k`1´tk x q " 0 λppxq|b| ppxq´2 b y`2 γpb y´∇y l k`1´tk y q " 0 (10) where ∇l " p∇ x l, ∇ y lq. If b x and b y are not zero, then: Substituting Equation (11) into Equation (10): sgnpb y qTb ppxq´1 y`2 γpb y´∇y l k`1´tk y q " 0 where T " λppxqpp . Note that: sgnpb y q " sgnp∇ y l k`1`tk y q So Equation (12) can be expressed as: sgnp∇ y l k`1`tk y qTb ppxq´1 y`2 γpb y´∇y l k`1´tk y q " 0 Unfortunately, we cannot obtain the explicit solution of the Equation (15). We can use the Newton method to get an approximate solution. If b y is solved, b x can be easily determined using Equations (11) and (13). The process is shown as Algorithm 1.

End End
Another problem is that in practice we don't know d in p(x). We have tested two ways to approximate d. One way is to use edge preserving filter (e.g., bilateral filter) to give an approximation of d and keep the exponent fixed during the iteration; Another way is to replace d with Gpl k`1 q during the iteration and Gp¨q represents the Gauss convolution operator, In most cases, both methods can generate similar prominent results. However, in some cases, dynamic approximation would give better results than fixed approximation because dynamic approximation can give a more accurate approximation of d along with the iteration. To illustrate this, consider the associated heat flow to problem Equation (3): where l TT and l NN are the second derivatives of l in the tangent and normal direction to the isophote lines respectively. From Equation (16), we have two critical conclusions: 1.
The illumination image, l k`1 , becomes increasingly smooth over time.

2.
Diffusion speed in the tangent direction is always faster than that in the normal direction.
The first conclusion conforms to the smooth assumption of the illumination image. If the illumination image has non smooth parts, then the second conclusion guarantees that the solution can preserve these parts. Thus, l k`1 continuously gets closer to d with calculation. However, the convergence proof of the algorithm is difficult since the exponent is changing during the iterations. If the exponent is fixed, the convergence proof can be directly obtained because the objective function is fixed and the iteration of split Bregman is monotone decreasing in the function values. The strict proof can be found in references [34,35]. If the exponent changes during the iterations, then the objective function changes as well. The convergence proof in references [34,35] cannot be applied here. However, we have tested numerous experiments and our algorithm did converges in all the tests. We leave it for further work.

Update:
2.4.4. Update l k`1 : which corresponds to the constraint L ě I ą 0. The process is shown as Algorithm 2.
(2)Given l k`1 and t k , update b k`1 by using algorithm 1.

Relation to Previous Methods
Let us revisit the model in Section 2.2. If we set ppxq " 2 and remove the constraint l ě i, our model is equivalent to homomorphic filtering [36]. Retaining the constraint l ě i and fixing p(x) = 2, it is similar to a random walk, Ng's model [17] and McCann algorithm [12]. Thus, our proposed model generalizes previous models.

Numerical Results
We present numerical results to demonstrate the efficiency of the proposed model and algorithm. For color images examples, we use HSV Retinex. We compare our proposed model with three state of the art methods, HoTVL1 [18], Ng's method [17] and multiscale Retinex [37].
For all the tests, the recovered reflectance of our model is: where L = exp(l) is the illumination function obtained from Section 2.4, and I = exp(i) is the original image. Note that the reflectance image is sometimes over enhanced, and we add the Gamma correction illumination to the reflectance image after decomposition. The Gamma correction of L with parameter s is: where W is the value of the white pixel. Parameter s was set to Section 2.2 in the tests. Thus: and the global framework of our proposed method is illustrated in Figure 6. correction illumination to the reflectance image after decomposition. The Gamma correction of L with parameter s is: where Wis the value of the white pixel. Parameterswas set to Section 2.2 in the tests. Thus: and the global framework of our proposed method is illustrated in Figure 6.

Synthetic Images
In this subsection, we set λ = 80, γ = 10 3 and w = 10 9 . Simulated illumination is added to the original texture images, as shown in Figure 7, with the numerical results shown in Figure 8. The recovered image following our proposed method is visually superior. We use signal to noise ratio (SNR) to measure the similarities between the original and recovered images, as shown in Figure 9. SNR from our proposed method is significantly superior to the other methods. We further use structural similarity index (SSIM) and CIEDE2000 color difference to measure the texture similarities and perceptual difference between the original and recovered images respectively, as shown in Tables 1 and 2. We can see from tables that our proposed method is superior to the other methods. We note that HoTVL1 failed in these tests. The main reason is that the assumption in HoTVL1 is piecewise constant and piecewise linear, which means that the shadow should be piecewise linear. However, this is not the case of these tests. The shadow part is almost a constant and also has sharp edges. Hence the result of HoTVL1 is not satisfactory.

Synthetic Images
In this subsection, we set λ = 80, γ = 10 3 and w = 10 9 . Simulated illumination is added to the original texture images, as shown in Figure 7, with the numerical results shown in Figure 8. The recovered image following our proposed method is visually superior. We use signal to noise ratio (SNR) to measure the similarities between the original and recovered images, as shown in Figure 9. SNR from our proposed method is significantly superior to the other methods. We further use structural similarity index (SSIM) and CIEDE2000 color difference to measure the texture similarities and perceptual difference between the original and recovered images respectively, as shown in Tables 1  and 2. We can see from tables that our proposed method is superior to the other methods. We note that HoTVL1 failed in these tests. The main reason is that the assumption in HoTVL1 is piecewise constant and piecewise linear, which means that the shadow should be piecewise linear. However, this is not the case of these tests. The shadow part is almost a constant and also has sharp edges. Hence the result of HoTVL1 is not satisfactory. structural similarity index (SSIM) and CIEDE2000 color difference to measure the texture similarities and perceptual difference between the original and recovered images respectively, as shown in Tables 1 and 2. We can see from tables that our proposed method is superior to the other methods. We note that HoTVL1 failed in these tests. The main reason is that the assumption in HoTVL1 is piecewise constant and piecewise linear, which means that the shadow should be piecewise linear. However, this is not the case of these tests. The shadow part is almost a constant and also has sharp edges. Hence the result of HoTVL1 is not satisfactory.

Natural Images
For all tests, we set λ = 80, γ = 10 3 and w = 10 3 . We begin with Andelson's checkerboard shadow image, as shown in Figure 10a. Region A looks darker than region B, although they have the same values. Figure 11 shows the reconstructed illumination and reflection images using Ng's, HoTVL1, multiscale Retinex and our proposed model. HoTVL1 and the proposed method produce superior results to Ng's method and multiscale Retinex. The recovered illumination using our proposed method contains less reflectance information than HoTVL1, e.g., the outline of the cylinder. Our proposed method also contains less shadow information in the reflectance image than other methods. Table 3 compares the recovered intensity values of the two regions for the four methods. The contrast of the marked areas using our proposed method is superior to the other three methods.  Table 3 compares the recovered intensity values of the two regions for the four methods. The contrast of the marked areas using our proposed method is superior to the other three methods.  Consider the degraded image shown in Figure 10b. Figure 12 shows the reconstruction for the  Table 3 compares the recovered intensity values of the two regions for the four methods. The contrast of the marked areas using our proposed method is superior to the other three methods.    Consider the degraded image shown in Figure 10b. Figure 12 shows the reconstruction for the four methods. Note that in this example, we adopted a Gamma correction step, as discussed above. Our proposed method has superior visual outcome. Ng's method suffers halo artifacts, e.g., near the edges of the tower and the roof of the building, which rarely appear in HoTVL1, multiscale Retinex and our proposed method. However, many fine structures lost in the HoTVL1 and multiscale Retinex reproduced image, which are retained in our proposed method.
The next illustrative example is recovery of non-uniform degraded images. The two images in Figures 10c and 3b suffer from the strong shadow areas. Figure 13 shows the comparison for the considered methods. The shadow is almost entirely removed by our proposed method, whereas the other methods retain partly shadowed regions.
In the end, we test the effect of different approximations of d. We use bilateral filter to approximate d and keep the exponent fixed during the iterations. Figure 14 shows the numerical results. We see in Figure 14 that the illumination and the reflectance of the results are not as good as those in Figure 11. This experiment supports our discussion in Section 2.4.2. (c) (d)

Conclusions
We proposed a variable exponent functional model for Retinex, proved the existence of the solution for the model and provided the theoretical derivation. The proposed method can be applied to general degraded cases as well. Experimental results validatethat our proposed method can remove non-uniform illumination and significantly reduce halo artifacts.