IMAGE RETINEX BASED ON THE NONCONVEX TV-TYPE REGULARIZATION

. Retinex theory is introduced to show how the human visual system perceives the color and the illumination eﬀect such as Retinex illusions, medical image intensity inhomogeneity and color shadow eﬀect etc.. Many researchers have studied this ill-posed problem based on the framework of the variation energy functional for decades. However, to the best of our knowledge, the existing models via the sparsity of the image based on the nonconvex (cid:96) p -quasinorm were limited. To deal with this problem, this paper considers a TV p -HOTV q -based retinex model with p,q ∈ (0 , 1). Specially, the TV p term based on the total variation(TV) regularization can describe the reﬂectance eﬃciently, which has the piecewise constant structure. The HOTV q term based on the high order total variation(HOTV) regularization can penalize the smooth structure called the illumination. Since the proposed model is non-convex, non-smooth and non-Lipschitz, we employ the iteratively reweighed (cid:96) 1 (IRL1) algorithm to solve it. We also discuss some properties of our proposed model and algo-rithm. Experimental experiments on the simulated and real images illustrate the eﬀectiveness and the robustness of our proposed model both visually and quantitatively by compared with some related state-of-the-art variational models.

1. Introduction. Our vision tends to see the same color in a given image regardless of the light and the color of objects remains relatively constant under varying illumination. The illusion effect is called the retinex effect. The word "Retinex" is made up of two words: retina and cortex. The retinex theory was initially proposed by Land and Mccann [27] through modeling the color perception of the human visual system, which can ascertain the reflectance of a object in which both reflectance and illumination are unknown. In general, we would like to remove the retinex effect to strengthen the image with the nonuniform lighting by separating the illumination from the reflectance in the given image. However, this problem is a classic ill-posed problem due to the lacking of some prior information and the uniqueness.
Many implementations and improvements of the retinex problem have been studied in the literature to overcome the ill-posedness. A path-based algorithm was put forward by Land and Mccann in [27], in which they considered the reflectance at each pixel depending on the multiplication of the ratios along random walks. There were other algorithms based on this approach [36,46,8]. The major difference is the path geometry among these path-based algorithms. Piecewise linear paths were used in [28], while double spirals and Brownian paths were used in [8]. But these methods need to tune a mass of parameters and are generally slow in performance. Recursive algorithms were proposed by using a recursive matrix calculation to replace the path computation in [29,38]. These algorithms are more efficient than a path-based approach. However, it is not clear how many iterations are required in the recursive method and the final results are strongly affected by the number of the iteration setting. Jobson et al. in [22,23] described the homomorphic filtering scheme for the retinex problem, where they modeled the reflectance as a low-pass filter. Meanwhile, a kernel-based retinex (KBR) [4] was proposed through using the expectation value of a weighted random variable with a suitable kernel function. This scheme can provide an analysis of the action of the KBR on contrast, produce a two-sided contrast modification, and enhance both under-and overexposed picture. Recently, Rizzi [3] and Palma et al. [1] presented other perception based on the color correction in the framework of variational techniques, which allow more flexible local control of the contrast adjustment and attachment to the data.
Recently, the variation-partial differential equation (PDE) based models also got a lot of attention in the image retinex problem, where the models depend on the suitable energy functional space. For the PDE-based models, researchers usually decompose the intensity of the original image into a product of the reflectance and the illumination, where the reflectance and the illumination are usually assumed to be the spatially smoothing and the piecewise constant respectively. With the help of these assumptions, Morel et al. in [40] solved a Poisson equation and Horn in [20] used the Laplacian to the sum s = l + r, where s, l, and r are respectively the initial image, illumination, and reflectance in the logarithmic domain, and applied threshold function to clip the peaks of |∆s| in order to get a Poisson equation. There are some similar approaches to the work in [5,39,40]. However, the numerical methods for these models severely depend on the Courant-Friedrichs-Lewy (CFL) condition when solving the corresponding PDE.
Different to above PDE-based methods, a lot of attention was paid to the variationbased models recently due to the flexible numerical methods and the good mathematical properties based on the choosing functional spaces. Based on the assumption that the illumination is smooth, Kimmel et al. [25] proposed a variational model for explaining how the human vision system perceives colors. However, their model did not consider the reflectance piecewise constant assumption. As one improvement, Ma and Osher [37] coupled the total variation (TV) with the nonlocal TV as the regularization term to the retinex problem and then proposed the Bregman iteration scheme to solve it. But this model did not consider the relationship between the illumination and the reflectance. It was highly hard to present some theoretical analysis such as the existence and the uniqueness of the solution. In [42], Ng and Wang established a TV model based on the prior of the reflection function, where some constraints are added to make the model more appropriate and reasonable. Liang and Zhang [34] further proposed the higher order TV-based L1 model (called the HOTVL1) for suppressing the staircase artifact generated by the TV regularization [37,42]. Zosso et al in [53] established a unifying model based on the nonlocal formulation. In addition, Wang and He [50] also proposed a variational model with barrier functionals in order to satisfy some constraints. They solved their model by an alternating minimization and the steepest decent method. By using the L 0 quasi-norm regularization to penalize the sparsity based on some transformations, Duan et al. [12] proposed to combine the Mumford-Shah model with the L 0 quasi-norm for the retinex problem. Gu et al in [18] recently proposed a detail preserving variational model for simultaneously estimating the illumination and the reflectance from an observed image without the log-transform. However, most of the existing variational models perform a log-transform on the original image in their models to turn the product form into the addition form. This operation can reduce the computational complexity and simulate the human eye perception of light intensity. Furthermore, these existing models due to the lacking of the efficient penalty for the noise prior are not robust when the image contains some noises. In order to overcome these drawbacks, Liu et al [32] established a new TV-based model (called the ETV), where the data-fitting term uses the exponent transform to replace the log-transform and one of the regularization terms is used to describe the noise prior and their proposed model can enhance the image details. Aiming to extend this work, this paper employs the non-convex p -quasi-norm for improving the ETV and the main contributions are summarized as follows.
• We propose a new image retinex model based on the nonconvex total-variation (TV)-type regularization, where the data fitting term and one of regularization terms based on the exponent-transform are used to describe image details efficiently, one regularization term based on the nonconvex total variation pquasinorm to penalize the sparsity of the piecewise constants and one regularization term based on the nonconvex high-order total variation q -quasinorm to penalize the piecewise smoothing regions. • Since the proposed model is nonconvex, non-smooth and non-Lipschitz, we employ the alternating minimization method, where the iteratively reweighted 1 algorithm and the alternating direction method of multipliers are used in the sub-minimization problems. We also give some theoretical analysises of the proposed model and the numerical algorithm. • Numerical comparisons demonstrate that our proposed model achieves competitive results for the image retinex problem with the different types.
The rest of this paper is organized as follows. In section 2, we introduce the iteratively reweighed 1 algorithm. In section 3, we present the proposed model and some theoretical analysis about the existence of the solution. The efficient numerical methods are proposed to solve the proposed model in this section. Section 4 presents some numerical results to demonstrate the effectiveness of the proposed method compared with some state-of-the-art retinex methods. Finally, we conclude our paper in Section 5.
Notations. For simplification, we set a gray image f ∈ X. Here X =: R M ×N is equipped with the inner product and the norm as . Here the forward difference operator D + specifically denotes as Similarly, we can define the backward difference operator as Using above symbols, we get the discrete version of the Hessian operator as ). In addition, ∀ x := (x 1 , x 2 ) ∈ Y and ∀ y := (y 1 , y 2 ) ∈ Y, we define their inner product x, y = x 1 , y 1 + x 2 , y 2 .
Then the 1 norm and the 2 -norm are defined in the space Y as For the simplification, we set · 2,2 as · 2 in the following. As an extension, the p -quasinorm is defined by for p ∈ (0, 1). Let Z := Y × Y and ∀ ρ = (ρ 11 , ρ 12 ; ρ 21 , ρ 22 ) ∈ Z, the q -quasinorm based on the space Z is defined by for q ∈ (0, 1). Here F denotes the Frobenius norm. Let div and div 2 to denote the adjoint operators of ∇ and ∇ 2 respectively, based on the divergence theorem −divx, u X = x, ∇u Y and div 2 ρ, u X = ρ, 2. Iteratively reweighted l 1 algorithm. This section recalls some numerical methods related to our proposed model and the numerical algorithm.
2.1. The 2 − p proximal problem. During last decades, how to solve the 2 − p proximal problem, as one subproblem used in the field of compressed sensing, pattern recognition and image restoration, has attracted great attention [7,9,54,33,31]. Formally, this problem can be written as where p ∈ (0, 1) and X depends on some spaces such as Y and Z. That is to say, we set X := 2 if x ∈ X and set X := F if y ∈ Y. Since the numerical method based on a different space in X is similar, we in the following consider the case of x ∈ X as To solve the problem (2), there are some great challenges for the theoretical analysis and the numerical implementation due to its non-convexity, non-smoothing and non-Lipschitz. Many schemes have been proposed to overcome these drawbacks such as the smoothing approximation method [9], the general iterative shrinkagethresholding algorithm [54,31], the active set method [15] and the iterative reweighed 1 minimization algorithm (IRLA) [7,26,35]. Here we mainly use the latter to solve the problem (2). The IRLA was first proposed by Candes et al. [7] to slove the problem (2). This algorithm consists of solving a sequence of the weighted 1 -minimization problems, where the weighted parameter used for the next iteration are computed from the value of the current solution. Specifically, solving the original problem (2) is approximately replaced by It is obvious that this problem is the classical 2 − 1 problem and the closed-form solution can be obtained by using the soft shrinkage scheme In the following, we summarize the IRLA to solve the problem (2) as Algorithm 1.
If setting ρ = 1 in Algorithm 1, that is to say, we use the fixed smooth parameter η in the approximated problem (3), the convergence analysis has been discussed in [10]. However, the generated sequence {x k } only converges to the solution of the smooth problem (3) rather than the solution of the original problem (2). In order to give the convergence analysis of Algorithm 1, we need to show the existence of the solution for the problem (2).
In fact, the objective function in the problem (2) is bounded from below and coercive, it is easy to deduce the existence of the solution in the problem (2) as done in [15]. Now we consider the convergence of Algorithm 1 related to the framework in [24].
Theorem 2.2. The sequence {x n } generated by Algorithm 1 converges to one solution of the problem (2).
Proof. Since the problem (2) is separable, we only consider the case at the pixel point (i, j). Based on the problem (3), we can deduce that Omitting the iteration superscript in the problem (3), we can rewrite G (x i,j , x i,j , η) as Set x i,j := x k+1 i,j in the norm · 2 and η := η k in above equation, it is easy to get Combining it with the inequation (4), we have In addition, since η k+1 ≤ η, we deduce . and then obtain that This result illustrates that the sequence G x k i,j , x k i,j , η k is non-increasing. Based on the bounded blow of the problem (3) and η k → 0 as k → ∞, we can conclude that the sequence G x k i,j , x k i,j , η k has a accumulation as lim Now we only want to show thatx i,j is one solution of the problem (3). In fact, we plug the point (x i,j ,x i,j , 0) in the problem (4), it is easy to find that the problem (4) is equivalent to the problem (3) for i = 1, · · · · · · M and j = 1, · · · · · · N .
Recently, Ochs et al. [43] extended it to solve various kinds of iteratively reweighed algorithms for the computer vision problems.
3. Total variation model. In this section, we first introduce our proposed model to the retinex problem and then give an efficient numerical method to solve it.
3.1. The proposed model. The intensity in the observed image is determined by two parts: the intensity of the illumination received at this point and the ratio of the illumination reflected from this point [21]. Assuming that the illumination of one object is inhomogeneous, it is natural to model the observation through inducing the varying illumination [30]. The varying illumination may have many different sources such as the varying distance based on a point lighting source, the uneven thickness of clouds via filtering sun lights, or the biased magnetic field generated in an MRI machine. Regardless of the sources, we furthermore assume that the illumination should vary smoothly. For simplification, this paper studies and discusses our proposed model based on a single channel-like gray-level image. Specially, the observation image f can be modeled by: where denotes the entrywise multiplication and η denotes the additive noise. Here we need two assumptions, which include that the source reflectance R ∈ (0, 1) is the piecewise constant and the illumination L ∈ (0, ∞) is the spatially smoothing.
In the model (6), the multiplication operation leads to many numerical difficulties to obtain R and L. One efficient method is of decoupling with them by using the logarithm transformation as v = r + l, where r = log(R) and l = log(L). With this manipulation, the model (6) can be rewritten as The problem (7) is the classic additive model and many schemes were proposed to deal with it based on some prior assumptions. Here we assume that η is the white Gaussian noise and then consider the following retinex model 1 Since the sequence G x k i,j , x k i,j , η k has a accumulation, we only set x i,j =x i,j to get the value of this accumulation for the functional G(·).
where α, β and γ are the positive regularization parameter, µ is the penalty parameter. The term f − e v 2 2 is used for the fidelity, ∇e v 1 , ∇r p 2,p and ∇ 2 l q FF,q are the regularization term, and v − r − l 2 2 is the penalty term. The last term here is to ensure the well-posedness of the model, where θ can be taken to be very small.
We remark the difference between the ETV used in [32] and the proposed model (8). It is obvious that the model (8) is the ETV if setting p = q = 1. However, as is well known that the image structure will be sparse based on some transformations such as the gradient-based operators and the wavelet-based operators. So it is more reasonable for restricting parameters p, q ∈ (0, 1) in order to preserve the sparsity when proposing the reconstruction model as done in (8). This fact can be observed that our proposed model (8) can give better results in the quality of the solution compared with the ETV from the numerical comparisons.
Proof. Since F(v, r, l) is bounded from below, the existence can be followed from its continuity and coercivity. Let w : is coercive and then the minimizing sequence w k for the functional F w k must be bounded. Boundary of the sequence in the space X implies the existence of a convergence subsequence of w kt . Using the continuity of the function F(w) in the closed image domain, above assertion is held.

3.2.
The alternating minimization algorithm. The problem (8) is a multivariable optimization problem, where the variables v, r, l are coupled together. In order to decouple with them, one efficient method is of using the alternating minimization method to transforming it into some easy and solvable subproblems. In general, the alternating minimization (AM) method, which is a rather old and fundamental algorithm, is attractive due to its simplicity and efficiency. This method is widely used in the linearly constrained optimization problems such as the machine learning, the image and signal processing, as well as other fields [17]. Specially, we can summarize this method as Algorithm 2 in the following.
Based on Algorithm 2, we can find that the minimization sequence {F (v n , r n , l n )} is monotonous and nonincreasing. Then we can deduce that there exists a subsequence of {F (v n , r n , l n )} to converge to an accumulation point of the problem (8).
In the following, we consider how to solve every subproblem in Algorithm 2.
3.2.1. Subproblem (9). This subproblem is nonconvex and nonsmooth due to the existence of the nonlinear function e v and the 1 -norm. One efficient method is of using the alternating direction method of multipliers (ADMM) through introducing some auxiliary variables to break the original problem into smaller pieces, each of which are then easier to be handled [17,2,13,11,6,16,47,51]. If the original problem is convex, the ADMM is closely related to many other methods such as the dual decomposition, the method of multipliers, Douglas-Rachford splitting, Spingarn's method of partial inverses, Dykstra's alternating projections; see for the work in [17]. To this end, we need to first introduce two auxiliary variables u, b and then transform the problem (9) into an equivalent form as Based on the framework of the augmented Lagrangian method, we need to rewrite (12) as the equivalent unconstrained optimization problem where Λ 1 and Λ 2 are the Lagrange multiplier, r 1 and r 2 are the positive penalty parameter. Using the ADMM, we can solve the saddle point problem (13) by the following Gaussian-Seidel iteration scheme In the iteration scheme (14)-(18), it is not easy to solve the subproblem (14) since the nonlinear function e v makes the objective function be nonconvex. So we either use an approximate solution as a substitute in the update which might cause divergence, or solve the subproblems by numerical algorithms which can bring computational burden. One efficient method is of using the Taylor linearizations around and then adds a proximal term as This is a smooth and convex optimization problem, the optimal condition of which can be obtained by Furthermore, we have To the subproblem (15), it is a convex and smooth. With the help of the variational method and simultaneously choosing the Neumann boundary condition, we can get the corresponding Euler-Lagrangian equation as follows (20) ( where I is the identity operator and ∆ is the Laplace operator. This problem is of solving the linear equation, and the numerical method depends on the boundary condition [52]. This paper mainly uses the periodic boundary condition. So the matrix operator (1 + r 1 )I − r 2 ∆ is a circulant matrix and the solution can be obtained by using the fast Fourier transform (FFT): where F denotes the Fourier transform and F −1 denotes its inverse transform.
To solve the subproblem (16), it corresponds to the classical 2 − 1 problem and the closed-form solution can be obtained by using the soft-shrinkage operator as (10). It is not easy to solve this problem due to its nonconvexity, nonsmoothing and non-Lipschitz. In order to overcome these numerical difficulties, we need to use the operator splitting scheme by introducing an auxiliary variable w := ∇r and then obtain a constrained optimization problem

Subproblem
where v n+1 := v k+1 . Here v k+1 denotes the outputting of the inner iteration related to the subproblem (9). Using the augmented Lagrangian method, we can obtain a saddle point problem as where Λ 3 is a Lagrange multiplier, and r 3 is a positive penalty parameter. Based on the ADMM, we can solve the problem (22) through the following iteration scheme Now we consider to solve the smooth subproblem (23), the optimal condition of which can be written as Based on the assumption of the periodic boundary condition again, the solution of this problem can be obtained by using the FFT.
To the problem (24), we rewrite it as

Subproblem (11).
The method for solving this problem is similar to solving the problem (10). So we similarly introduce an auxiliary variable z := ∇ 2 l and use the ADMM to solve it as follows.
In the iteration scheme (29)-(31), the subproblem (29) is smooth and convex. Then the optimization condition can be written as (µ + θ)I + r 4 div 2 ∇ 2 l = µ v n+1 − r n+1 + r 4 div 2 z ω + div 2 Λ ω 4 . Using the assumption of the periodic boundary condition, the solution can be obtained by using the FFT through Similarly, the subproblem (30) can be written as (33) min which is the 2 − q problem. Then we still use Algorithm 1 to solve it by setting in the problem (1).

Remark 1.
To Algorithm 2, there include some inner iterations to solve the subproblem (9)- (11). Specifically, we use the ADMM to solve them respectively. In general, more exact solution can be obtained when increasing the number of iteration. However, we notice that a better decomposition result can be still obtained by setting the inner iteration to be one. So we choose it as the stopping condition for all of the inner iterations in the numerical implementation when solving our proposed model.

Numerical experiments.
To demonstrate the performance of our proposed model (8), this section presents a series of numerical implementations to remove the different illusions and noisy levels of the different kind of the images. To this end, we choose some state-of-the art models such as the TVH1 [42], the HOTVL1 [34], the L0MS [12] and the ETV [32] compared with our proposed model (8).
As is well known, choosing suitable parameter in the proposed model is essential for obtaining a good correction and the reconstruction result. For the numerical algorithm, the suitable parameters can improve the convergence rate and the numerical stability. However, there involve more than two parameters such as the regularization parameters α, β and γ, the penalty parameters µ, p and q. In general, three ways can be proposed for how to select them. One is of trying on some values for providing satisfactory results and then keeping it to be the fixed outputting. The second way is by using Game theory to reformulate the optimization problem so that the dimension of the parameters' space is much reduced, as is done in [48]. The third approach is of choosing them by using some rules based on the prior statistical characteristic such as the L-curve method [19], the generalized cross validation (GCV) [14,49], the discrepancy principle [41,45], or the variational Bayes approach [44]. In our experiments, we choose the former approach. Specifically, we set α ∈ [0.01, 0.04] according to the severity of noises, β ∈ [0.001, 0.02] depending on the illusion, µ ∈ [0.7, 1.2] and p ∈ [0.9, 0.97] based on the image structure. In addition, we fix γ = 0.08, τ = 1e − 5, q = 0.85 in the numerical implementation. In addition, in order to quantitate the performance, we adopt the Peak Signal to Noise Ratio (PSNR) and the Mean Structural Similarity (MSSIM) to evaluate the effectiveness of the used models. All comparison algorithms are implemented and tested in Matalb R2018a under Windows OS on a laptop computer equipped with Intel Core i7 2.7GHz. 4.1. Synthetic images. We first start with reconstructing two synthetic images degraded by the white Gaussian noise with σ = 0.001 and a suitable bias field as shown in Figure 1. From the values of the PSNR and the MSSIM in Table 1, we can conclude that our method demonstrates competitive score. Since humans are involved in evaluating the quality of the reconstructed result, it absolutely becomes essential to provide visual quality of the reconstructed images. As shown in Figure 2, we can observe that the recovered reflectance and illumination images generated by our proposed model are less noise effects and artifacts than other models. Especially, some noises can be found in Figure 2 (a), (b) and (c) and some piecewise constants can be found in Figure 2    In the medical imaging, the obtained images may be corrupted by the bias field due to non-uniform illuminations or the parallel magnetic resonance imaging (MRI) 3 . The correction of the bias field is similar to the retinex problem where we need to remove the light effect caused by the illumination. Here, we use two slices of T 1 -weighted brain MR image as examples, which are contaminated the white Gaussian noise with different standard deviation as 0.03, 0.05, 0.07, 0.09 and the intensity non-uniformity as 40%. The quantitative comparisons as shown in Table 2 indicate that our proposed model achieves the best results in terms of the PSNR and the MSSIM.
To show the effectiveness via the visual quality of the retinex images, Figure 3 shows the related comparisons between our proposed model and other models. It is easy to observe that all models can provide visual preferable results compared with the original images. Specially, removing the noise by using our proposed model and the ETV is more robust than other models. In order to show the robustness of our proposed model, we zoom a part of the retinex image degraded by the noise with the standard deviation as 0.09 as an example and then plot them in the bottom column of Figure 3. Again, we still observe that the ETV and our proposed model can reconstruct more reasonable images. The important difference is that our proposed model gives more natural structures such edges and the piecewise constant regions.

Intensity inhomogenous images.
Here we select one slice from the T 1weighted brain volume and then generate ten degraded images by adding the white Gaussian noise with the standard deviation σ = 0.001 and the different intensity inhomogeneities. The PSNR and the MSSIM for the reconstructed images by using different methods can be found in Figure 4. It obviously concludes that our proposed model outperforms other four models based on these quantitative criterions. Meanwhile, it can be observed that our proposed model is more stable than other models. That is to say the change of the intensity inhomogeneities can not seriously effect the reconstruction effectiveness of our proposed model due to the gentle change of the PSNR and the MSSIM as shown in Figure 4. In general, the coefficient of variation (CV) is a standardized measure of the dispersion of the probability or frequency distribution. Then the CV can be used to estimate the effectiveness of the different method for reconstructing the degraded image. Formally, the CV of each tissue T defined by  Table 2. PSNR and MSSIM of T 1 -weighted brain MRIs with different levels of Gaussian white noises.
where σ(T) and µ(T) are the standard deviation and the mean of the intensities in the tissue T. Here we employ it for evaluating the White Matter (WM), the Gray Matter (GM) and the Cerebrospinal Fluid (CSF) on the bias corrected images.
The related boxplots are shown in Figure 5. It is similarly observed that the CV values based on our proposed model outperform other models, which implies that our proposed model can effectively deal with noises and intensity inhomogeneities again. Now we consider to compare the visualization through choosing three degraded images as shown in Figure 6. The reconstructed images and the different images between the clean images and the restored images are shown in Figure 7. We can find that all methods can remove the noises and different intensity inhomogeneities efficiently. Especially, we can observe the colorbars of the different images based on our proposed method more dark and less structural information compared with other models, which implies the effectiveness of our proposed method. In addition, since the relative error between the clean images f and the restored images f 1 can also display the effectiveness of the different model, we still use three images in Figure 6 as examples and then plot these curves of R(f, f 1 ) as shown in Figure  8. It is obvious that R(f, f 1 ) based on our proposed model is smaller than other methods. These illustrate the advantage and effectiveness of our proposed model again.
4.1.3. Convergence analysis. Since our proposed model (8) is nonconvex and non-Lipschitz, it is not easy to prove the convergence of the proposed algorithm. In order to show the convergence from the numerical implementation, we need to visually analyze the relative error of the solution r n and l n R(r n , l n ) = r n − r n−1 2 2 r n 2 2 , l n − l n−1 2 2 l n 2 2 and the energy of the objective function in the problem (8). In general, R(r n , l n ) tending to zero implies the convergence of the numerical algorithm. Here we only use the third image as shown in Figure 6 as an example. As we can see from Figure 9 (a) Input that R(r n , l n ) and the numerical energy decay when increasing the iteration, which demonstrates that our proposed algorithm converges well numerically.

Real MRIs.
In this experiment, we present to reconstruct some real MRIs such as the bladder and brains as shown in Figure 10, and the heads as shown in Figure 13, where these real MRIs are obviously corrupted by the noise and the bias field. The reconstructed results presented in Figure 11 demonstrate that our proposed model can well recover the real intensities and improve the quality of the original MRIs. In addition, in order to show the reconstruction effectiveness, we select a slice of each image which locations are shown in Figure 10, and then   Figure 7. Performances of five methods on T 1 -weighted brain images with different intensity inhomogeneities for the first, third and fifth rows; Colorbar to different between the clean images and restored images for the second, fourth and sixth rows.
plot their gray value curves as shown in Figure 12. Based on Figure 12 Numerical Energy (c) Figure 9. The relative errors of r and l and numerical energy of our model for the third image in Figure 6.
can observe that the TVH1 and the HOTVL1 can not remove the bias field very well and the L0MS model makes the image too smooth. In contrast, our proposed model and the ETV can remove bias field and noises effectively. Furthermore, by the careful observation, the distribution of the gray values based on our proposed model is relatively closer to the original image than the ETV. That is to say, the reconstructed results generated by the ETV reduce the contrast of the images. In other words, our proposed model outperforms other models. Meanwhile, we can get similar conclusion from Figure 12(b) and (c). Compared with the MRIs as shown in Figure 10, the MRIs of the head as shown in Figure 13 have more structure information. By observing the contour images of the original MRIs from the bottom of Figure 13, it is not easy to find the detailed information due to some effects of the bias field and the noise. These effects can reduce the effectiveness of the subsequent image processing such as the image segmentation and the image identification. Hence, it is very important to correct the bias field and the noise. In the following, we consider to reconstruct these MRIs. The reconstructed images and their contour images are respectively shown in Figure  14 and Figure 15. Figure 14 shows that the reconstructed comparisons generated by the different models can present visually preferable results compared with the original images. To intuitively observe the reconstructed results, we display the contours of the reconstructed results as shown in Figure 15. It is easy to observe some facts as follows. The TVH1 is able to remove non-uniform illuminations effectively, but it can break the background. The L0MS cannot recover the bottom of the corrupted images well. Compared with the HOTVL1 and the ETV, our proposed model not only removes the bias field more robust, but also preserves the details of the images.  Figure 16. For the Adelson's checkerboard shadow image, though visually region B is brighter than region A, they are of the same intensity value. From the comparative results in Table 3, it can be observed that all of the PSNR and the MSSIM of the proposed model are higher than other models. In addition, the reflectance component R and the illumination component L obtained by different models are displayed in Figure 17. The visual results demonstrate that our proposed model can well decompose the reflectance and the illumination.

5.
Conclusion. In this paper, we presented an efficient model to decompose the reflectance and the illumination based on the TV p -HOTV q regularization. Since the proposed model is non-convex, non-smoothing and non-Lipschitz, we employ the iteratively reweighed 1 algorithm based on the alternating minimization algorithm to solve it. Some mathematics properties of the proposed model and the numerical method were discussed. Various numerical results were implemented to demonstrate the advantages of our proposed model over some state-of-the-art retinex methods.