A TV Bregman iterative model of Retinex theory ∗

A feature of the human visual system (HVS) is color constancy, namely, the ability to determine the color under varying illumination conditions. Retinex theory, formulated by Edwin H. Land, aimed to simulate and explain how the HVS perceives color. In this paper, we establish a total variation (TV) and nonlocal TV regularized model of Retinex theory that can be solved by a fast computational approach based on Bregman iteration. We demonstrate the performance of our method by numerical results.


Introduction
The image obtained by a digital camera is different from the one perceived by the human visual system (HVS). Digital images record information based on the light reflected by objects while human being can automatically discount the variation of the illumination. This feature is called color constancy which ensures that the perceived color remains constant under varying illumination conditions. Edwin H. Land's Retinex theory [13,14,15,16] is the first computational model that aims to simulate the HVS. The word "Retinex" is a portmanteau of retina and cortex because the eye and the brain are involved in the system. The basic assumption of Retinex theory can be summarized as follows: (1) The HVS performs the same computation in each of three independent color channels (RGB).
(2) In each color channel, the image intensity I is propotional to the reflectance R and the ilumination E: I ∼ RE (3) The reflectance R can be perceived by the HVS while the illumination E is automatically discounted by the HVS.
In practice, in each color channel, we suppose that the image intensity at pixel x is the product of the reflectance and the illumination.
To simulate the HVS, we need to recover the reflectance R from a given image I.
For example, in Figure-1, the image is called "Adelson's checker shadow illusion". Area A seems darker than area B for most people, but actually they are the same as shown in the right image, i.e., I(A) = I(B). We notice that area B is in the shadow of the cylinder, so the illumination of area A is larger than that of area B, i.e., E(A) > E(B). By Retinex theory, the reflectance of area A is smaller than that of area B, i.e., R(A) < R(B). For human being, the illumination is discounted by the HVS. So area A seems darker than area B. The recovering of reflectance is a mathematically ill-posed problem and there are many formulations for Retinex theory. In this paper, we introduce the total variation (TV) and nonlocal TV regularizers which make the original problem well-posed. These problems can be solved by an efficient approach based on the Bregman iteration. This paper is organized as follows: In section 2, we briefly review several different implementations of Retinex theory. In section 3 we present the TV and nonlocal TV regularized model and computational algorithms based on the Bregman iteration. In section 4, numerical results are shown to demonstrate the algorithm's performance. In section 5, we discuss the future work for the proposed model.

Previous works
Retinex algorithms are basically categorized [17,12,19] as path-based algorithms, recursive algorithms, center/surround algorithms, PDE-based algorithms, and variational algorithms. We will briefly review these algorithms in this section.

Path-based algorithms
Land's original works [13,15] consider the reflectance at each pixel depending on the multiplication of the ratios along random walks. Brainard et al. [2] described the path algorithm with stochastic theory. These kinds of methods need a large number of parameters and have high computation complexity.

Recursive algorithm
Frankle and McCann [5,6] replaced the path computation by a recursive matrix calculation. It highly improved the computational efficiency. But the number of iterations is not clearly defined and can strongly influence the final result.

Center/Surround alogrithms
The center/surround approach was proposed by Land, et al. [13] and later improved by Jobson, et al. [11]. Jobson et al. then proposed the SSR (single scale Retinex) and the MSR (multi-scale Retinex). The idea is that the output reflectance values can be computed by substracting a blurred version of the input image. This algorithm is easy to implement but needs many parameters.

PDE formulation
In PDE-based formulations, people usually take a logarithm for the original formulation to get i = r + e where i = log(I), r = log(R), and e = log(E). A basic motivation is that the reflectance corresponds to the sharp details in the image (i.e. edges) whereas the illumination is spatially smooth. Horn [10] applied the Laplacian to get ∆i = ∆r + ∆e and then use a threshold function to clip the peaks to get a Poisson equation ∆r = τ (∆i). The threshold function is defined by Blake's [1] and Morel's [18,19] formulation are similar to Horn's.
For example, Morel took a gradient first to get ∇i = ∇r + ∇e. By the assumption that ∇e is relatively small, he applied the threshold function componentwise to ∇e and then took the divergence to get Poisson equations.
where ǫ depends on the size of the image.
The Poisson equation can be solved by fast algorithms such as by an FFT, and there is only one parameter t in the threshold function τ .

Variational formulation
By the same assumption in PDE formulations, Kimmel et al. [12] present a variational Retinex formulation.
subject to: u ≥ i and ∇u, − → n = 0 on ∂Ω. This is a quadratic programming problem that can be solved by many methods such project normalized steepest descent method as in Kimmel's paper. However, this regularizer is not as effective as total variation in the image restoration problems.

Proposed algorithm
In this section, we present a total variation regularized formulation and an effective computational algorithm based on the split Bregman method. We discuss the connection of our formulation with the previous PDE formulations. We also introduce a nonlocal TV regularized formulation.

A TV regularized model
The total variation regularizer introduced by Rudin, Osher, and Fatemi [22] is very effective in recovering edges of images. This coincides with the basic motivation as in PDE-based algorithms: the reflectance corresponds to the sharp details in the image (i.e. edges) and the illumination is spatially smooth. So we can formulate the problem as a TV regularized minimization problem.
Our formulation can be written in the following form: The first term is the TV regularizer term that ensure us to find the sharp details. The second term is the L 2 term of the gradient of the illumination. The motivation of this term is the smoothness of the illumination. The parameter t balances the two terms.
We can adapt this model to gamma-corrected images. The gamma-correction applies a concave function to the raw images, in practice a logarithm or a s γ power with 0 < γ < 1. Assuming that the gamma-correction is logarithmic (s γ has a similar shape as logarithm), we can replace i = log(I) by I and have the gamma-correction model

Bregman iteration
Bregman methods [3] are effective for solving sparse reconstruction problems in image processing [20]. Bregman iteration is based on the definition of Bregman distance. The Bregman distance for a convex function J is given as: where J is convex but not necessarily differentiable, such as the L 1 norm, and H is convex and differentiable with zero as its minimum value. This problem can be solved by Bregman iterations: Equation (3.4) is usually referred as "adding back noise". The advantage of Bregman iteration is to transform a constrained problem into a series of unconstrained subproblems. And the unconstrained subproblems (3.3) can be solved by the split Bregman method.
The split Bregman was introduced by Goldstein and Osher [8] for solving L 1 , TV, and related regularized problems and applied to various imaging problems [9,23]. The split Bregman method aims to solve the unconstrained problem: where J and H are as before, and Φ is linear operator. For instance, in ( The key idea of the split Bregman method is to introduce an auxiliary variable d = Φu, and try to solve the constrained problem For simplicity, we introduce a new variable b k = p k d /λ. We notice that p k d = λb k and p k u = −λΦ T b k , and thus the iterations become: d k+1 and u k+1 can be updated alternatively. We first fix u k to update d k+1 and then fix d k+1 to update u k+1 . Then we get the general split Bregman iterations.

The TV Bregman iterative algorithm
We can apply the split Bregman method to our variational problem by setting Φu = ∇u, J(Φu) = t Ω |∇u|, and H(u) = 1 2 ||∇(u − i)|| 2 2 . By the split Bregman technique, we introduce auxiliary variable d = (d x , d y ) such that d x = ∇ x u and d y = ∇ y u and get the constrained problem: By the general split Bregman method (3.5, 3.6, 3.7), we have The subproblems (3.8) and (3.9) can be explicitly solved [8]. We summarize them in the following algorithm.
The isotropic shrinkage function is defined by In iteration (2), we update u k+1 by solving a Poisson equation associated with the zero Neumann boundary condition. The solution of this equation is not unique. But the difference of any two solutions is a constant. By fixing the value of one pixel, we can get a unique solution. In practice, the equation is solved by a discrete cosine transformation.

Connection to PDE-based algorithms
In the first iteration, we can solve for d 1 = shrink t (∇i). We need to find u 1 by solving the Poisson equation. ∆u = divd 1 = div shrink t (∇i) This is very similar to the PDE-based algorithms mentioned above. The only difference is the threshold function. In previous formulations, the threshold function clips the peaks and thus produces some jumps. Here we shrink the peaks so that the result is still continuous.

A nonlocal TV regularized model
Following Buades, Coll, and Morel's nonlocal mean method [4], Gilboa and Osher introduced the nonlocal TV regularizer [7]. This is another successful method in image processing, especially for textured images [23].
For an image u(x) : Ω → R, we can define the nonlocal weight between two pixel x and y: where G a is the Gaussian kernel with standard deviation a. With nonlocal weights, we can define the nonlocal gradient operator as the vector of all partial difference ∇ w u(x, ·) at x such that: And the nonlocal TV regularizer can be defined as So the nonlocal TV regularized model for Retinex theory is as follows: And the corresponding gamma-correction model is Similarly, we have a numerical algorithm using the split Bregman method. In the discrete case, we write images u and i as vectors by connecting the columns from left to right and suppose the size is n, i.e., u, i ∈ R n . Then the gradient operator can be represented as a 2n by n matrix, denoted by where D x = ∇ x and D y = ∇ y are two n by n matrices. For the nonlocal case, the nonlocal weights come from the original image i, or I if using gamma-correction model. In practice, we choose m closest patches for each pixel and set their weights 1 and weights of others 0 (see [7,23]). Then ∇ w can be written as an mn by n matrix: where D i is an n by n matrix with the diagonal elements -1 and one nonzero element 1 in each row. So d, b are vectors in R mn We apply the split Bregman method to the nonlocal TV regularized model by introducing an auxiliary variable d = ∇ w u. In discrete version, d can be written as d = (d 1 , d 2 , · · · , d m ) with d j = (d j,1 , d j,2 , · · · , d j,n ) ∈ R n . And the minimization problem is When we plug this into the general form of the split Bregman iterations, we have We explicitly solve the equations (3.11) and (3.12) and get the nonlocal TV Bregman iterative algorithm: Here shrink t is the nonlocal isotropic shrinkage function. For v ∈ R mn , denote v = (x 1 , x 2 , . . . , x m ) T and

Numerical results
In this section, we will show implementation details and some numerical results. Here we use the gamma-correction model for the following examples. We also implement the gammacorrection version of Morel's [19] PDE-based algorithm, i.e., replace i by I, and compare the results.
Before calculation, we perform color balance by linearly stretching the input image range to [0,255] for each color channel. The color of resulting images is balanced in the same way. The threshold parameter t depends on images. In each iteration, we calculate ||u k − u k+1 ||/||u k+1 || and stop iteration if it is less than tolerance ǫ = 0.02.
For the nonlocal TV regularized model, we first need to calculate the nonlocal weights from an image f and the nonlocal TV operator. For each pixel i = (i 1 , i 2 ), we choose pixel j = (j 1 , j 2 ) in a (2w + 1) × (2w + 1) window centered at i and calculate the weight w(i, j) = e −d(i,j) where d(i, j) is the distance between i and j: where G a is the Gaussian with standard deviation a, and p is called the patch size. After calculation, for each pixel, we choose m largest weights and set them 1 and set the others 0. Then we can obtain nonlocal TV operator D w as in the discussion in previous section. In our implementation, the nonlocal weights come from the input image I and we set w = 10, a = 3, p = 5, and m = 8. In the following examples, we test different t and choose the best one for each image and each algorithm.
The first example in Figure- (121, 103, 107). The contrast of the latter image is stronger than that of the former one. Since we use gamma-correction model, the difference between the original image and resulting image is the illumination by the assumption. We find that the difference image from the TV regularized model is just the shadow of the green cylinder while the difference image from the PDE-base model carries more information of the object itself. In the second example in Figure-3, the image is a piece of text in a shadow of some object. By the color constancy assumption, we are supposed to get a constant background after computation. We compare the results of PDE-based model, the TV regularized model, and the nonlocal TV regularized model. In the resulting image of the PDE-based model, the shadow is partially removed and the background is not constant. The shadow in the resulting image of TV regularized model is less obvious and the background is more constant than that of the PDE-based model. The resulting image of nonlocal TV regularized model is better than the other two images. We can barely see the shadow and the background is constant almost everywhere. The last example in Figure-4 is a piece of cloth with colorful bands and each band is of a constant color. In the image, the right part is under a shadow. After the computation, we are supposed to remove the shadow. In the resulting image of the PDE-based algorithm, the shadow is partially removed and the contrast of the difference image is not strong. The shadow in the resulting image of the TV regularized model is weaker than that of the PDEbased model, but it is not fully removed. If we choose a larger t, the shadow is weaker, but the difference between similar color is eliminated too. Actually, the difference between similar color is weakened. The result of the nonlocal TV regularized model has almost no shadow and each color band is constant. And its difference image has the strongest contrast among these three.

Future work
The TV and nonlocal TV regularized model can be modified to the problem of image enhancement. Our approach can be simply generalized to 3D to process 3D medical image such as MRI data. We can also apply our model to hyperspectral images.