A COUPLED TOTAL VARIATION MODEL WITH CURVATURE DRIVEN FOR IMAGE COLORIZATION

,


1.
Introduction. Image colorization is the process to turn a grayscale image to a color image, which is always a challenging subject in the fields of image processing and computer vision. It has been widely used in the film industry to make old productions more attractive, and has been lately used in the medical imaging technique to enhance the visual recognition of medical images [27]. The reverse problem called decolorization also attracts much attention in recent years [7,12]. In the literature, there are mainly two categories of colorization methods, namely, color propagation and color transfer.
In the color propagation approach, the colorization process is carried out by propagating the input color data to the whole grayscale image. The methods for color propagation include optimization, diffusion, variational techniques, see for instance [8,11,9,29,25]. Although this approach is efficient and often generates visual appealing results for colorization, the initial color of grayscale image need to be input manually. Another colorization approach [28,18,3,23] is based on color transfer. In this approach, the colorization process is performed by transferring colors from a source color image to a target grayscale image. Currently, some exemplar-based colorization techniques [28,3] are used for color transfer. These techniques usually make the assumption that pixels with similar intensities or similar patches should have similar colors. Generally, the exemplar-based methods suffer from spatial color consistency problems since the colorization of each pixel is processed independently. In addition, the selection of source color images is very important. Recently, a unified approach [24] by combining both exemplar and input color is proposed for image colorization.
In this paper, we focus on the colorization method based on color propagation. Let us briefly discuss existing color propagation methods. Levin et al. [9] considered the assumption that neighboring pixels with similar intensities should have similar colors. Based on this idea, they proposed to solve an optimization problem to diffuse the chrominance of input color pixels to the whole grayscale image. Yatziv and Sapiro [25] proposed to inpaint color pixel values by minimizing the difference between the gradient of luminance channel Y and the gradient of chrominance channel Cb (or Cr) in the YCbCr color space. Subsequently, Sapiro et al. [29] developed a fast image colorization approach by using geodesic distances to blend the chrominance of input color to the other grey-value pixels. The above-mentioned methods do not consider the preservation of contours of grayscale images. In order to preserve image contours in the coloriztion process, some variational models in luminance channel (or brightness channel) have been proposed for colorization. Kang et al. [8] first proposed a coupled variational model for image colorization. This method can preserve the contours quite well. Since this variational model is not convex, its performance can be affected by initial guesses of color images significantly. Kang et al. [11] further developed a fast image colorization method using vector-valued reproducing kernel Hilbert spaces. However, the model is still non-convex even the efficiency can be enhanced. Some other coupled total variation methods for colorization can be found in [4,10]. It is interesting to note in existing coupled variational methods that color diffusion is controlled by using the gradient information of target grayscale image. Even the gradient driven diffusion can preserve the contours of target grayscale image quite well, the pseudo-edges appearing in the contours can incorrectly affect the color diffusion process. This is the main disadvantage when we would like to colorize a large region of grayscale image.
In this paper, we propose a new coupled total variation model based on curvature information (or isophoto curve information) of target grayscale image to control the colorization process. Moreover, we formulate the optimization model for colorization in the YCbCr color space. There are two distinct advantages of this approach. The proposed model makes use of curvature information to control the color diffusion process more effectively for colorizing some large regions in grayscale images. On the other hand, the involved optimization problem is convex and several convex programming solvers can be employed to determine the solution which is not sensitive to initial guesses of color images.
The rest of this paper is organized as follows. In Section 2, we present and study the proposed method. In Section 3, the existence of the minimizer of the proposed model is shown, and the algorithm is given to solve the resulting model. In Section 4, experimental results are reported to show that the performance of the proposed model is competitive with the other color propagation methods for testing images. Finally, some concluding remarks will be given in Section 5.
2. Variational methods. Throughout this paper, let Ω ⊂ R 2 be the bounded image domain with Lipschitz boundary. Write D ⊂ Ω as the domain where we wish to colorize, and D c = Ω \ D where the color pixels are given. In the RGB color space, denote by a target grayscale image to be mapped with some colors given in D c . We remark that when (x, y) ∈ D, i.e., no color is given. For a grayscale image f : (x, y) ∈ Ω → R, we denote by the gradient and the Laplacian of f respectively. It is known that f is usually used to describe curvature information of the isophoto lines in f [26]. There are some typical color spaces such as RGB space, two color-difference components YUV (luminance and chrominance) space, YCbCr (luminance and chrominance) space and CB (chromaticity and brightness) space, see for instance [17,19]. Most of colorization systems work in luminance-chrominance spaces such as YUV space or YCbCr space. In this paper, we present our method in the YCbCr color space.

2.1.
A non-convex model. In [8], Kang et al. presented a variational method for large region colorization via weighted harmonic map. They gave the method in the CB color-space which is a non-convex color space. Notice that the luminance component of target image, B 0 = |f 0 | is known everywhere in Ω and the chromaticity component C 0 = f 0 /B 0 are only known in D c , where |·| denotes the Euclidean norm of vectors. If the chromaticity data C in D can be gotten form C 0 by some color propagation scheme, a final RGB image f will be recovered by f = CB 0 . Therefore, they proposed the following variational functional to recover the chromaticity C and the corresponding minimization problem where g 0 : R + → R + is a monotone decreasing function such that g 0 (0) = 1, g 0 (t) > 0 for any t > 0, and lim t→+∞ g 0 (t) = 0. The function g 0 in [8] is taken as g 0 (t) = e (−0.002t 2 ) . Here is the Gaussian kernel with A, σ > 0. In (2), S 2 is given as follow: {Z ∈ R 3 | |Z| = 1}. Thus the model in (2) includes a non-convex constraint |C| = 1. Hence, Kang with α > 0 and the corresponding minimization problem (4) min{F α (C) : C ∈ W 1,2 (Ω; R 3 ), |C| ≤ 1 for a.e. x ∈ Ω}.
Furthermore, they proved that the solutions {C α } α of the variation problem (4) converge to a solution C of the variation problem in (2) in terms of Γ-convergence as α tends to infinity. They also employed gradient descent method to numerically solve the approximation model in (4) when α is large enough.
There is a disadvantage in the variational method by gradient driven diffusion for large region colorization. Let us discuss it by considering an example on the colorization result of Kang's model (2). In Figure 1, image (a) denotes the target grayscale image with given color and image (c) is the plot of diffusion function g 0 (|∇B 0 |). We see from the original image (b) that the desired colorization result should be that the left pepper of image (a) is fully colorized by red color. It is clear from image (c) that in middle of the left pepper there exists an interior pseudo-edge where the values of diffusion function g 0 are close to 0. This leads to that the given color cannot further diffuse to the whole region of left pepper and stops at this pseudo-edge. It follows that the left pepper is locally colorized by Kang's model as shown in image (d). On the other hand, another disadvantage for Kang's model is that the model is not convex. This results in that the colorization results of Kang's model depend on the initial guesses (see Figures 5 and 6), and some current fast algorithms based on convex optimization cannot been directly applied to solve the model.

2.2.
The proposed method. In this paper, we first present our method in a convex color space: YCbCr color space to ensure that the proposed objective function for colorization is convex. Our colorization method can be explained as follows. For a target grayscale image f 0 ∈ L ∞ (Ω; R 3 ) with some colors given in D c , the luminance data Y 0 ∈ L ∞ (Ω) is known everywhere in Ω and (Cb 0 , Cr 0 ) ∈ L ∞ (D c ; R 3 ) are known in D c by using the transform map from RGB color-space into YCbCr color-space. If we can get the chrominance channels (Cb, Cr) values in Ω form (Cb 0 , Cr 0 ) by some color propagation scheme, a final RGB image will be recovered by combining the luminance channel Y 0 . Therefore, the colorization problem can be reformulated as the problem of reconstructing the two chrominance channels. On the other hand, we would like to overcome the effect of some pseudo-edges in target image which hamper the given color to propagate further. Here we introduce a new color driven diffusion and propose a coupled total variation model for colorizaztion such that a large region colorization is achieved possibly. More precisely, the following variational problem is studied to recover the chrominance channels: and g is a monotone decreasing function, is the laplacian operator as defined in (1), G σ is the Gaussian kernel as described in (3) and λ is a fidelity paramter to balance the two terms in the above objective function. After the minimizer V is determined, the color image can be obtained by combining the luminance Y 0 and the recovered chrominance channels V. Here we conclude the proposed method for colorization as follows: 1) For an image f 0 with color given in D c , compute Y 0 and V 0 by the transformation form RGB space to YCbCr space.
2) Recover V by using the Algorithm 1 designed later to solve the proposed model (5). 3) The reconstructed image f is computed from (Y 0 ,V) by the inverse transformation form RGB space to YCbCr space.
As an example, we show in Figure 1 the colorization results of the proposed method. We first present in image (f) the colorization result of the proposed model in which g(| G σ * Y 0 |) is replaced by g(|∇G σ * Y 0 |), and image (e) is the plot of the diffusion function g(|∇G σ * Y 0 |). It is clear from image (e) that there is similarly an interior pseudo-edge as described as in image (c). This leads to that this model also cannot colorize the whole region of left pepper; (see image (f)). Then, we give the colorization result of the proposed model in image (h), and image (g) is the plot of the diffusion function g(| G σ * Y 0 |). We see form image (g) that the interior pseudo-edge described by laplace controlling is not so clear as that in images (c)-(e) by gradient controlling. It follows that the given color can cross the pseudo-edge to further propagate into the whole region of left pepper, and then stops at the external boundary. We can see this fact from image (h). These results indicate that the proposed model can effectively overcome the effect of some pseudo-edges in target image such that large region colorization may be achieved.
3.1. Theoretical analysis. In the proposed model (5), the controlling function g is chosen as With this choosing, we can give the definition of g-BV(Ω; R 2 ). In fact, regarding the weighted bounded variation (WBV) space, the detailed discussions can be found in [1,2,6]. It is direct to extend the space to the vectorial case. Here we only present some basic results on the lower semi-continuity and compactness for the theoretical analysis purpose.
Let Ω be a bounded open subset of R n , and g(x) ≥ 0 be a continuous real valued function on Ω. Define g-BV(Ω; R s ) (s ≥ 1) as a space of functions V ∈ L 1 (Ω; R s ) such that the following quantity Lemma 3.3. Assume {V k } is a sequence in g-BV(Ω; R s ) satisfying sup k u k g-BV < ∞. Then there exists a subsequence {V kj } ∞ j=1 and a function V ∈ g-BV(Ω; R s ) such that V kj → V in L 1 (Ω; R s ), when j → ∞.
Next we show the existence of the solutions for (5). Let Ω ⊂ R 2 be a bounded domain with Lipschitz boundary. From now on we always write the energy functional to be E(V) in (5). The theoretical analysis of the proposed model is given as follows. Proof. It is clear that E(V) ≥ 0 for all V ∈ Λ. In addition, we can construct a functionV such thatV ∈ g-BV(Ω; R 2 ) and E(V) < +∞. Therefore, we consider a minimizing sequence {V n } ∞ n=1 ⊂ Λ in (6).
By the definition of the minimizing sequence, we get that there exists a constant M > 0 such that E(V n ) ≤ M . That is, On the other hand, it follows form V n ∈ Λ that V n L ∞ (Ω;R 2 ) ≤ √ 2b for all n ∈ N. Thus, we conclude that {V n } is uniformly bounded in g-BV(Ω; R 2 ). By using the compactness of g-BV(Ω; R 2 ), we deduce that, up to a subsequence, {V n } converges to V * = (Cb * , Cr * ) ∈ L 1 (Ω; R 2 ) and does so a.e. in Ω, i.e., As a consequence of lower semi-continuity of g-BV(Ω; R 2 ), Recalling that V n L ∞ (Ω;R 2 ) ≤ √ 2b for all n ∈ N, by using Lebesgue dominated convergence theorem and (8), we obtain Hence, we conclude Meanwhile, it is obvious that a ≤ Cb * (x), Cr * (x) ≤ b for a.e. x ∈ Ω. Thus, V * is indeed an optimal solution of the variational problem (5). (5) is convex. The uniqueness of solution refers that there may be several global optimal solutions attaining the same minimum objective function value.

Remark 1. The objective function in
3.2. The convex programming solver. In this subsection, we present the basic procedure to solve the proposed model (5). We remark that this model is convex in the setting of YCbCr color space. Our idea is to use the alternating minimization by incorporating the variable-splitting technique [15]. By introducing the notations Γ ≡ {V = (Cb, Cr)|a ≤ Cb, Cr ≤ b}, and the auxiliary variables d ≡ (d 1 , d 2 ) and W, we reformulate (5) as the following constrained optimization problem: subject to ∇V = d, V = W and W ∈ Γ.
The above optimization problem is a standard convex program on convex set with linear constraints. Therefore, it fits the framework of the alternating direction method of multipliers (ADMM), which is introduced in [20,14,22]. Attaching the Lagrangian multiplier λ = [ λ 1 T , λ 2 T ] T to the linear constraints, the augmented Lagrangian function of (9) is given by where µ 1 , µ 2 > 0 are penalty parameters. χ Γ (W) is the indicator function given by The final algorithm of problem (9) is presented as follows: Algorithm 1.

Go back to step 2 until
Notice that the convergence of ADMM [13,5] can be used here for the proposed Algorithm 1 for solving the problem in (9). We conclude the convergence of Algorithm 1 as the following theorem.
where f true and f are the true color image and the colorized image, respectively. When the PSNR of the colorized image is high or its MSE is low, it refers to a good colorized image generated by the method. In the following tests, Y 0 is employed and computed via the finite difference method instead of G σ * Y 0 . The convolution of G σ and Y 0 is used for the purpose of mathematical analysis. The parameter τ in the controlling function g (7) is set to be 1000. The fidelity parameter λ of the proposed model is chosen by providing the highest PSNR value of the colorized image. In Algorithm 1, we set the stopping criterion to be = 0.5 × 10 −4 , and the penalty parameters are fixed and set as follows: µ 1 = 15, µ 2 = 0.75. Indeed, we find that the numerical performance of ADMM is not very sensitive to the choices of µ 1 and µ 2 , provided that they are not extremely large or small. All the computation is done under MATLAB implementation on a personal computer with an Intel(R) Core(TM) i3-2370M CPU and a 2.40-GHz central processing unit. , Yatziv and Sapiro [29], and Kang and March [8].
The source codes of the Levin's method and Yatziv's method are from the websites 1 2 , respectively. In the following tests, we use their default settings to generate colorized images. For Kang's method, the gradient descent method is used to solve the model (4), the stopping criterion is C k+1 − C k ≤ 0.001. The parameters for Kang's model (4) are chosen as follows: α is set to be 10 as described as in the paper [8]. The fidelity parameter λ is selected corresponding to the highest PSNR value. In [8], the authors pointed out that the use of λ = 10 is good enough for the method.  Figure 2.
In the first experiment, we show the colorization results of the three state-of-theart methods and the proposed method. In Figure 2, it is clear from images in Figure  2(c) and 2(d) that the color of type I is given (24.07%) more than that of type II (3.47%). It is seen from the images in Figure 2(h) and 2(l) that the recovered colors by Yatizv's method have been changed in the colorization process. When the area of the given color is large (the image in Figure 2(e)), the colorized images by Levin's method, Kang's method and the proposed method are very good, see the images in Figure 2(g), 2(i) and 2(j). However, when the area of the given color is small (the image in Figure 2(f)), the colorized images by Levin's method, Yatizv's method and Kang's method are not good, see the images in Figures 2(k), 2(l) and 2(m). Both Levin's method and Yatizv's method cannot preserve the contours of the target image. Although Kang's method can preserve the contours well, it is obvious from the image in Figure 2(m) that there are some areas which cannot be colorized. The colorized image generated by the proposed method is very good, see the image in Figure 2(n). In Table 1, we show that the PSNR value and the MSE value of the colorized images by the proposed method are quite good compared with the images generated by the other methods, in particular when the given color area is small.

4.2.
The effect of parameters. In this subsection, we study several parameter settings in the proposed method, namely, the fidelity parameter, the penalty parameters and the initial guess of the algorithm. Here, we use the images in Figures 2(e) and 2(f) to test the effect of these parameters.
We first investigate the performance of the proposed model with respect to the fidelity parameter λ. We observe form Figure 3 that the PSNR value is monotonically increasing when λ varies form 1 to 6, and the PSNR value is slightly decreasing after λ arrives at 10. These results indicate that the performance of the proposed method is slightly affected for large values of λ. In Kang's model, we see from Figure 4 that the performance of Kang's method is about the same for a wide range of λ.
On the other hand, we find that the proposed method is not sensitive to initial guess, whereas Kang's method can be sensitive to initial guess. In Figure 5 and 6, we show three different initial guesses and their effect to the resulting colorized images. When the given color region is large, the performance of Kang's method is very good. However, when the given color region is small, its performance is dependent on different initial guesses. In Tables 2 and 3, we further show that the PSNR and MSE values of the colorized images by Kang's method are sensitive to different initial guesses. In Figure 7 and 8, we further demonstrate that this sensitivity cannot be ignored even the parameter λ is allowed to tune. We see from the figures that for a given initial guess, the PSNR values of the colorized images for different values of λ by Kang's method are about the same. These results may show the disadvantage of using a nonconvex model for colorization. In contrast, we show in Figure 5

Experiment 2.
In this subsection, we show the performance of different methods by colorizing some testing grayscale images in [8,9], see the first row of Figure 9. The colorization results by different methods are given in Figure 10 where the given color regions in grayscale images are shown in the second row of Figure 9. The percentages of the number of color pixel values are set to be between 8% and 11%. For Yatizv's and Levin's methods, we use their default settings to generate colorized images. For Kang's and the proposed methods, we set the parameters that can give the largest PSNRs of the colorized images.
In Table 4, we report the PSNR and MSE values of the colorized images in Figure 10. We note that the performance of the proposed method is competitive with those of the other methods. We see from Figure 10 that the colors of all the colorized images by Yatizv's method are changed mildly. The performance of Kang's method is not good especially for images I and II. Also the colorized images by Levin method is not good compared with those by the proposed method, for example, see image I. For the colorized images I, II and III, the colorization performance of the proposed method is quite good because the total variation model with curvature driven is particularly useful to colorize large regions with clear contours.
In Figure 11, we display some zoomed parts of colorized images by different methods and compare them with the original color images. We find that the generated colors by the proposed methods are close to the original colors. For example, let us look at the zoomed region of an "orange" in Image I. We see that the colors of Levin's, Yatizv's and Kang's methods are not matched with the original colors. However, the colors generated by the proposed method are very close to those in Kang's method The proposed method initial guess (0,0,0) (1,1,1) (1,0,0) Table 4. The values of PSNR and MSE of the colorized images in Figure 10, and the corresponding value of λ used in the proposed method.
color Image I. In Image II, we can compare the colorization results in the zoomed region of a "flower". It is clear that there is a sharp boundary in the original color image. The colorized images by Levin's and Yatizv's method are not good. The red color of the flower also appears in the other regions, i.e., the mechanism of colorization is not controlled well. In contrast, both Kang's and the proposed methods can capture the boundary of the flower and give the red color to the flower suitably. However, Kang's method fails to colorize appropriately in the outside region of the flower. We see that the performance of the proposed method can detect boundary of the flower and colorize both inside and outside regions with suitable colors. In Image III, we see that the proposed method is very good in colorization of red colors to a "pepper". By comparison with Images I, II, III and IV, we see that Image IV does not have clear boundaries inside the picture. Therefore, we do not expect that the performance of the proposed method is superior to the other testing methods. Finally, in Table 5 we show the computational times of different methods. We remark that the program for Yatizv's method is to use C++ version, while the others are to use Matlab implementation. The computational time of Yatizv's method is less than those of the other methods. According to the tables, we see the the proposed method is quite competitive with the other methods in terms of colorization and computational performance.   5. Concluding remarks. In this paper, we have proposed a coupled total variation model based curvature driven for image colorization. The existence of solutions of the proposed model have been shown in the paper. We employed an effective algorithm to solve the variational model numerically, and we showed the convergence of this algorithm. Experiments results have been given to illustrate the effectiveness of the proposed method, and its performance is quite well compared with the other testing methods.
Image I Image II Image III Image IV Figure 9. The first row: original color images; the second row: corresponding grayscale images with their percentages of given color regions to be 9.76%, 10.34%, 10.60% and 8.10%.
Original Image I: Original Image II: Original Image III: Figure 11. The zoomed regions of the colorized images by different methods. Each first row is the original color image. Each second row contains the colorized images: Levin's method (1st column); Yatziv's (2nd column); Kang's method (3rd column) and the proposed method (4th column).