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.

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

Figure 1

Elementary structure of the retina, lateral geniculate nucleus, and cortex.

Figure 2

Receptive field in the retina. (a) on-center ganglion cell; (b) off-center ganglion cell.

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 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:

The reflecting object a Lambertian reflector and reflectance corresponds to sharp details in the image;

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.

2. New Assumption and Proposed Model

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

Figure 3

Images with different illumination and the corresponding surfaces: (a) text image; (b) book image; (c) surface of (a); (d) surface of (b).

Figure 4

A single row extracted from: (a) text image; (b) book image.

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.

2.2. 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→∞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 , 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:

Figure 5

Schematicdiagram for Retinex.

We first convert Equation (2) into the logarithmic domain: so that:

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 , 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, near the edges of d where the gradient is large, and so the regularizer is similar to a TV regularizer which can preserve edges; 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.

2.3. Solution Existence

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

(variable exponent spaces): Let be a bounded open set with Lipchitz boundary and be a measurable function, with the family of all measurable functions on being . We define a functional, which is also called modular: and a norm:

Then the variable exponent Lebesgue and Sobolev spaces are, respectively: and:

With the norm , becomes a Banach space.

(relationship between modular and norm [ be a modular on and , then +1.

(embedding theorem [ , and for a.e. . Then is continuously embedded in .

(convexity [ , with as in model (3). Then for each , is convex in .

(weak lower semi continuity [ be bounded from below, and the map is convex in each . Then the energy functional, , is weak lower semi-continuous in .

Let be a bounded open set with Lipchitz boundary, , then the minimization problem: has a minimizer .

Let , be the minimizing sequence for . Then: where M denotes a universal positive constant that may differ from line to line. Hence . Thanks to Poincare inequality, we have , and from Lemma 2, . Therefore, , and together with the inequality , we obtain . This implies that is a uniformly bounded sequence in due to Lemma 1, and is also uniformly bounded in . Since is a reflexive Banach space, up to a subsequence, there exists such that converges weakly to in . From Lemma 4, is lower semi continuous in . Thus:

Therefore, is the minimum point of . ☐

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

2.4.1. Subproblem l with Fixed b and t

Given the fixed variable b and t , our aim is to find the solution of the problem: which has the optimality condition: where and . 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.

2.4.2. Subproblem b with Fixed l and t

Similarly, we solve: which has the optimality condition: where . If and are not zero, then:

Substituting Equation (11) into Equation (10): where . Note that:

So Equation (12) can be expressed as:

Unfortunately, we cannot obtain the explicit solution of the Equation (15). We can use the Newton method to get an approximate solution. If is solved, can be easily determined using Equations (11) and (13). The process is shown as Algorithm 1.

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 during the iteration and 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 and are the second derivatives of in the tangent and normal direction to the isophote lines respectively. From Equation (16), we have two critical conclusions:

The illumination image, , becomes increasingly smooth over time.

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

2.4.4. Update :

which corresponds to the constraint .

The process is shown as Algorithm 2.

2.5. Relation to Previous Methods

Let us revisit the model in Section 2.2. If we set 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.

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

Figure 6

Global framework of the proposed method.

3.1. Synthetic Images

In this subsection, we set λ = 80, γ = 103 and w = 109. 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 Table 1 and Table 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.

Figure 7

Synthetic example. (a,d) original image; (b,e) simulated illumination; (c,f) synthetic image.

Figure 8

Image reconstruction by different methods. (a,e) Ng; (b,f) HoTVL1; (c,g) proposed; (d,h) mutiscale Retinex.

Figure 9

SNR curves. (a) For Figure 7a; (b) For Figure 7d.

Table 1

SSIM of the four methods.

Image	Ng	HoTVL1	Proposed Method	Multiscale Retinex
Figure 7a	0.9076	0.6342	0.9289	0.8764
Figure 7d	0.7603	0.6633	0.8291	0.7238

Table 2

CIEDE2000colordifference of the four methods.

Image	Ng	HoTVL1	Proposed Method	Multiscale Retinex
Figure 7a	26.8927	29.6729	26.6182	27.9038
Figure 7d	24.5474	26.0115	21.9972	25.4528

3.2. Natural Images

For all tests, we set λ = 80, γ = 103 and w = 103. 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.

Figure 10

Test images. (a) Checkerboard image; (b) Tower image; (c) Text image.

Figure 11

Illumination (a–d) and reflection (e–h) images. (a,e) Ng’s method; (b,f) HoTVL1; (c,g) proposed method; (d,h) multiscale Retinex.

Table 3

Recovered intensity values for regions A and B of Figure 10a.

Image		Original	Ng	HoTVL1	Proposed Method	Multiscale Retinex
Checkerboard	A	120	140	85	135	109
	B	120	180	174	230	149

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.

Figure 12

Reproduced degraded image (Figure 10b). (a) Ng’s model; (b) HoTVL1; (c) proposed model; (d) multiscale Retinex.

The next illustrative example is recovery of non-uniform degraded images. The two images in Figure 10c and Figure 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.

Figure 13

Reproduced shadowed images (Figure 10c and Figure 3b): (a,e) Ng’s model; (b,f) HoTVL1; (c,g) proposed model; (d,h) multiscale Retinex.

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.

Figure 14

Reproduced illumination and reflectance (Figure 10a) (a) Illumination with fixed approximation; (b) Reflectance with fixed approximation. Note the circled part by the red square.

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