WAVELET FRAME BASED COLOR IMAGE DEMOSAICING

. Color image demosaicing consists in recovering full resolution color information from color-ﬁlter-array (CFA) samples with 66.7% amount of missing data. Most of the existing color demosaicing methods [14, 25, 16, 2, 26] are based on interpolation from inter-channel correlation and local geometry, which are not robust to highly saturated color images with small geometric features. In this paper, we introduce wavelet frame based methods by using a sparse wavelet [8, 22, 9, 23] approximation of individual color channels and color diﬀerences that recovers both geometric features and color information. The proposed models can be eﬃciently solved by Bregmanized operator splitting algorithm [27]. Numerical simulations of two datasets: McM and Kodak PhotoCD, show that our method outperforms other existing methods in terms of PSNR and visual quality.

1. Introduction.Single sensor (CCD/CMOS) are often used in conventional digital cameras to sample three color components like red, green and blue.Through color filter array (CFA) the sensor records one color per pixel location.In order to get a full-resolution color image, it is essential to interpolate the two missing color components at each pixel location.Since the CFA sample system has mosaic pattern for different colors, this process of recovering full-resolution images has been widely called as "color image demosaicing".
The most widely used CFA pattern is so-called Bayer pattern [1] (see Figure 1).It measures green channel on a quincunx grid with subsampling rate 50% , and red and blue channels on rectangular grids with subsampling rate 25%, due to the fact that the human vision is more sensitive to the wavelengths around the green light.
Color image demosaicing can be essentially regarded as an image inpainting or interpolation problem.If we denote Λ r , Λ g and Λ b as the index sets for the observed red, green and blue pixels, respectively, f = (f r , f g , f b ) as the observed color vector where P is the usual projection operator.For simplicity, we rewrite the equation (1) in vector form as (2) P Λ u = f where P Λ is the projection onto the tensor index sets Λ = Λ r × Λ g × Λ b .
Recently, several wavelet frame based inpainting algorithms have been developed in literature (see [5,4,3,10,6]), where the advantages of sparse wavelet frame approximation based methods were shown in terms of sharp edges and fine geometric features.This motivates us to further explore the application of wavelet frame in color image demosaicing problem.
1.1.Existing methods.Essentially, demosaicing can be seen as an vector-valued image interpolation problem with many possible solutions.For recovering a natural and good quality color image, it is important to exploit the inter-channel correlation instead of interpolating three channels independently.In literature, one way to utilize the inter-channel correlation is based on the constant-hue assumption, which assumes that in local area, the ratio or difference of two color channels are relatively constant or smooth.This is illustrated in Figure 2 for an example image.In addition, we apply hard threshold on the wavelet coefficients of the example image, i.e., set those coefficients less than 10% of the range as 0, then plot the histogram in Figure 3.In fact, the red-green different image has 20% more zero coefficients (after hard threshold) compared to the red channel image in wavelet domain.Based on this observation, it is possible to interpolate the missing pixels in red and blue channels using relatively over-sampled green ones.For example, the well-known Hamilton-Adams [16] method first interpolates the green channel by taking into account the second order derivatives of the red/blue channel, and then recover the red and blue by interpolating the channel difference u r − u g and u b − u g .The interpolation is efficient for images with low color saturation and smooth chromatic gradient.Later on, many work were done to improve the accuracy of inter-channel difference estimation, such as [17,24,19].
In [25], interpolation error was treated as "demosaicing noise" and a filtering procedure was exploited to remove the noise.This idea was further improved in [21] where a spatial adaptive method was proposed to eliminate the demosaicing noise.Similar idea was adopted in [2,26], where a non-local patch based approach was applied for noise elimination.Besides interpolation in spatial domain, many approaches are based on frequency domain.One consideration of frequency-domain approach is the application of adhoc diamond-shape low-pass filter [12].Yet as we know that using low-pass filter cannot accurately keep sharp edges and tiny features.In recent years, researchers have proposed methods combining both spatial and frequency reconstruction [14,18].For example, the alternating projection method [14], is a combination of spatial domain interpolation and projection onto a constraint set in the wavelet domain.
1.2.Our approach.Most of existing demosaicing methods are based on either spatial smoothness or high color correlation, and they may fail when these conditions are not satisfied.In addition, as observed in several work [2,26], natural images often present abrupt color changes and relatively high saturation, which makes interpolation difficult without introducing any artifact.Therefore, we aim to propose a regularization type model which can simultaneously recover sharp edges and small important structures in each color channel while maintaining the color consistency.More precisely, we apply a sparse reconstruction model on both color components and channel differences using tight frame system [8,22,9] in two stages.The main advantage of wavelet frame based approaches is that piecewise smooth images can be sparsely approximated by properly designed wavelet frames, and hence, the 1 -norm regularization of frame coefficients can bring us sparse representation which are closed to real images.Tight frame systems have shown good performance in various image restoration tasks, such as [5,4,3,10].For this particular application, we design a two-stage model in order to gradually recover the fine features and color information.
There are two sparsity prior models in image restorations: analysis based approach and synthesis based approach.It is a well known fact that analysis based approach is suitable for relative smooth images while synthesis based approach recovers features and edges better (see e.g.[11]).We provide an easy criterion to automatically choose one of the two approaches based on mean color saturation.Essentially, based on a crude initial guess, we determine which method should be applied according to the Mean Saturation (MS) level.In particular, we test on two standard data sets : McMaster (McM) and Kodak PhotoCD.For Kodak data set , 21 out of 24 images have MS less than a given threshold, and the analysis approach outperforms the synthesis one.For McM images, 15 out of 18 have MS greater than the threshold, and synthesis approach performs better.This coincides with the fact mentioned above, since it is known that images in Kodak data set are relatively smooth, while McM images have more features and sharp edges.This also shows numerically that MS used in this paper is a simple and reasonable criterion in practice.As a result, our method can outperform the best available methods at 2.6dB at most and 0.5dB at least on average for the total 42 test images.Finally, although our proposed method is iterative, it is fast taking into account of the quality of reconstructed images, for example, the LDI-NAT method of [26] has the closest PSNR value to our method, however, it takes more than half an hour to process a 500×500 image, while ours takes less than one minute under the same environment.This paper is to introduce wavelet frame based method into the area of demosaicing and to show that it can be used to improve the quality of demosaicing.
The rest of paper is organized as follows.Section 2 gives a brief introduction of tight frame based sparse recovery models.The main method and algorithms for color image demosaicing will be described in Section 3. Finally, numerical results on two data sets will be presented in Section 4 in comparison with other methods along with Conclusions in Section 5.
2. Tight frame based sparse recovery models.
2.1.Sparse recovery models.It is well known that piecewise smooth images can be sparsely approximated by wavelet tight frames which can be numerically computed by 1 -norm regularization of frame coefficients.For the simplicity, the tight frame transform can be written in matrix form, forward transform matrix W and inverse transform matrix W T , the tight frame property implies that W T W = I, while the identity W W T = I does not necessarily hold.In computation, we do not use the matrix multiplication.Instead, we use the fast wavelet frame decomposition and reconstruction algorithms (see r.g., [9]), a brief introduction of wavelet tight frame will be given in subsection 2.2.
Next, we introduce two models: the analysis based approach and the synthesis based approach.The analysis based approach for tight frame demosaicing can be formulated as a constrained minimization problem (3) min where • 1 means 1 -norm sum on each component, µ is the weight parameter.Meanwhile, the synthesis based approach solves the following problem (4) min where d is the coefficients of u in the tight frame transform domain.
As pointed out before, the analysis based approach tends to recover smooth images with fewer artifacts, while the synthesis based approach tends to explore more sparsity in the transform domain, hence, it recovers features and sharp edges better.We might apply different models for different images.
2.2.Wavelet tight frame.We briefly introduce the concept of wavelet tight frame here.Interesting readers can consult [23] for a short survey on theory and applications of frames, and [11] for a more detailed note.
A countable set where If X is a tight frame, then X is called a wavelet tight frame and Ψ is called wavelet.The multiresolution analysis (MRA) based wavelet can be generated by the unitary extension principle (UEP) of [22].In particular, we use B-spline wavelet frame of orders (0, 1, and 3) constructed in [22].For example, the 0th order of B-spline wavelet frame is a constant function and the corresponding wavelet tight frame is nothing but Haar wavelet frame [15], the other two are piecewise linear and cubic spline wavelet frame.Given a 1-D wavelet frame system for L 2 (R), the corresponding 2-D wavelet tight frame system for L 2 (R s ) can be constructed via the tensor products of 1-D wavelet frame (see e.g.[8,11]).The discrete wavelet transform can be generated by the filters of wavelets and corresponding refinable function that generate the MRA.A discrete image u is an 2-D array.We will use W to denote fast tensor product wavelet frame decomposition operator and use W to denote the fast reconstruction operator.Then by the UEP [22], we have W W = I, i.e. u = W W u for any image u.When the multiple level decomposition is used, we will further denote an L-level wavelet frame decomposition of u as W u = {W l,i,j u : 1 ≤ l ≤ L, (i, j) ∈ I}, where I denotes the index set of all frame bands.More details on discrete algorithms of wavelet frame transforms can be found in [11].
The wavelet frame system is a redundant system.The redundancy can reduce the amplification of possible error such as noise and artifacts generated during reconstruction process.More details on discrete algorithms and their relative analysis of wavelet frame transforms can be found in [11].
3. A two-stage wavelet frame based method.In this section, we propose a two-stage color image demosaicing method.First of all, as used in several existing demosaicing methods, color differences are used to interpolate unknown color channels since they are smoother than individual channel in most of natural images.Instead of applying regularization independently on each channel u r , u g , u b , we first apply regularization on the green channel u g and the inter-channel differences u rg := u r − u g and u bg := u b − u g to take advantage of color correlation.According to saturation degree, we apply automatically either sparsity-promoting synthesis based approach or smoothness-promoting analysis based approach.More precisely, synthesis based approach is applied to get highly sparse approximation for images with high saturation and abrupt color changes, and analysis based approach is more suitable for those having low saturation with smooth chromatic gradient.At the second stage, we apply a finer regularization on each channel u r , u g , u b and the inter-channel differences u rg , u bg , u br := u b − u r for correcting and recovering small structures while avoiding false colors.Here is the details of the overall method: 1) Stage 1.For a given image, we apply one of either analysis based model or synthesis based model.Which one to use is according to the mean saturation value of the initial guess.The details on the criterion of choosing models by the mean saturation value of an initial guess will be given in Section 4.
Next, we describe the two specific models in details.For convenience, we redefine the vector u = (u rg , u g , u bg ) T .Let µ = (µ rg , µ g , µ bg ) T be a regularization weight vector.The analysis based approach model is define as: T , then the above model can be rewritten as (7) min Similarly we can define the synthesis based approach as where d = (d rg , d g , d bg ) are the coefficients corresponding to u rg , u g , u bg , respectively.The recovered color image is obtained by transforming the coefficients vector back to image space u = W d and then sum up to get u r , u g , u b .For both models, we observe that a low order wavelet filter gives a sufficient good result.The system we used is undecimal Haar wavelet frame, for which 1-D filter is After the first stage, the image quality can be significantly improved, but we can still observe color artifacts.In order to refine small features and improve color consistence, we add more constraints and apply synthesis based approach on each color channel and the inter-channel differences.Let d = (d r , d g , d b , d rg , d bg , d br ) T be the coefficients of u r , u g , u b , u rg , u bg , u br in tight frame transform domain respectively and µ = (µ r , µ g , µ b , µ rg , µ bg , µ br ) T be the weight parameter.We consider the following regularization model: Similarly, we can reformulate the above problem in a compact form as (10) min .
Note that we do not need to store the whole matrix of B, since the matrix vector multiplication and its adjoint can be easily implemented.In this step, in order to reconstruct finer structures, we apply piecewise linear undecimal wavelet frame with 3.1.Algorithms.To implement these two steps, we need to solve the synthesis based approach ( 8) and (10), and analysis based approach (7).For the past few years, many efficient algorithms have been developed for solving these two classes of problems arising in image processing and compressive sensing, such as Bregman based methods [20,13] and other augmented Lagrangian based, splitting based methods.Among them, we will apply Bregman operator splitting algorithm (BOS) proposed in [27] since it can maximally decouple the variables when A T A is not identity matrix.The BOS is a method based on Bregman iteration and forwardbackward operator splitting [7], which can transfer the constrained problem into several easy and efficient subproblems without inner iterations.Another advantage of the algorithms is that it can be easily implemented in parallel which is attractive for large size image processing problems.The BOS method can be also interpreted as an inexact Uzawa method under the primal dual framework [28] and it can be applied in a large variety of inverse problems.For more applications and convergence proofs on this type of method, one can refer to [27,28].
In the following, we will directly apply the BOS method on both synthesis based and analysis based approaches.We take the synthesis based model ( 8) as an example.Starting from an initial guess u 0 , we set It is well known now that the second subproblem can be efficiently solved by the pointwise "shrinkage" operation: The problem (10) is solved in an analogous way by replacing the matrix A by B and the variables accordingly.For analysis based approach (7), we apply the BOS method and obtain the following iteration scheme: Since W T W = I, the second subproblem has a close formula and the third one is also solved by the shrinkage formula (12).4. Numerical results.In this section, we show the numerical results performed on the simulated mosaic images from two standard test data sets: McM dataset introduced in [26] and Kodak PhotoCD1 .Due to space limit, only a subset of 8 images for each dataset is shown in Figure 4.
For the choice of parameters in our method, the two parameters δ and λ in the algorithms (11) and ( 13) are fixed as 1 according to the convergence requirements of algorithms.Like self-similarity driven demosaicing [2], our initial guess u 0 is obtained by Hamilton-Adams scheme.In the first stage, the 1 -norm of the frame coefficients d 0 rg 1 , d 0 g 1 , d 0 bg 1 are calculated to define the weight vector, µ := µ g ( for balancing different components.Finally, there is only one free parameter µ g to be tuned.On the other side, we note that theoretically the choice of µ g does not affect the results since the models (7), ( 8) and (10) remain the same.Therefore the choice of parameters for our proposed method is rather easy and stable in practical use.In Figure 7, we compare the proposed two-stage method with the sparse regularization on the three channels independently, i.e, we apply directly the analysis model ( 3) with u = (u r , u g , u b ).The figure shows that our proposed method can improve image quality and suffer less color artifacts due to the usage of inter-channel correlation.
In the following, we explain how to choose between two approaches in Stage 1 based on the mean saturation value of u 0 = (u 0 r , u 0 g , u 0 b ).At a given pixel u 0 i,j , we adopt the usual saturation definition: | (i,j) , then the average over all pixels is defined as the mean saturation (MS) s.Finally, given a threshold τ , we apply synthesis based approach if s ≥ τ , otherwise analysis based one.For Kodak image set, there are 21 out of 24 images with MS less than 0.4, while for McM set, 15 out of 18 images have MS greater than 0.4.Our numerical   results have show that the synthesis approach performs better on McM while the analysis approach performs better on Kodak.This demonstrates that MS level provides a simple criterion for automatically determining the sparse recovery model in practice.
The numerical simulations show that our method (TFD) is the only one which performs well on both McM and Kodak data sets.The average PSNR value of all 42 images from both datasets is given in Table 1, which shows that our method   2. Our method is much better than other methods and the only one close to ours is the LDI-NAT method, but it is much slower and performs poorly for the Kodak images.The zoom-in results are given in Figure 8, 9 and 10.Furthermore, Figure 8 shows that only our method can reconstruct the fine yellow line without huge color variation and in Figure 8 we observe almost no color artifacts in TFD recovered image.Finally, the PSNR values for the 8 Kodak images from Figure 4(b) are given in Table 3.The zoom-in results of one test image are also given in Figure 11.We can see that LPA outperforms TFD in terms of PSNR, while visual quality is similar.Moreover, the performance of LPA is much lower than TFD on McM images (almost 3dB difference).
Finally, in terms of computation time, our iterative methods are obviously slower than non-iterative ones, but it is much faster than patch based LDI-NAT method [26], for a McM image with 500×500 pixels, it takes our methods at most 40 seconds while more than half an hour for LDI-NAT method on a laptop with 2Ghz dualcore intel CPU and 2GB memory.It should be pointed out that the BOS method can be further accelerated by either applying an adaptive parameters strategy or implementation in parallel, which can potentially reduce the computation time in real application.As we mentioned early that our focus here is to introduce the wavelet frame based method into color image demosaicing, since the wavelet frame based method has been proven efficient (see e.g.[11]).parts and sharp at the edges while keep color consistence.Based on the mean saturation of an initial guess, our proposed method can automatically choose the better approach out of analysis based approach (7) or synthesis based approach (8) for the first step.In the second step, by enforcing the constraints and increasing the order of the wavelet filter, we can furtherly improve the image quality.It is shown in the experiments that the proposed approaches did generally better than the existing approaches on two different kinds of test data sets.The proposed method potentially can be applied for demosaicing problems with the presence of noise due to the nature of regularization.In future, we will explore more color statistics and regularization terms to automatically recognize and protect more types of features so that the color image demosaicing method can be improved furtherly.

Figure 1 :P
Figure 1: Bayer pattern used in single-sensor digital camera.

Figure 2 :Figure 3 :
Figure 2: A demo of inter-channel correlation of Boat image in 4(b).

Figure 5 :
Figure 5: Visual quality improvement from Stage 1 to Stage 2.

Figure 6 :
Figure 6: The improvement of PSNR over 8 test images from each data set (Figure 4).

Figure 7 :
Figure 7: Comparison of sparse regularization with and without inter-channel correlation.

Figure 8 :Figure 9 :
Figure 8: Zoom-in part of the demosaicing result for image 2 in McM image set.

Figure 10 :Figure 11 :
Figure 10: Zoom-in part of the demosaicing result for image 8 in McM image set.

Table 1 :
total average PSNR value (dB) of all the 42 images of two datasets.
5.Conclusions.This paper presents wavelet tight frame based method for color image demosaicing problem.By applying regularization on both color channels and inter channel differences, we allow the reconstruction result smooth in the cartoon

Table 2 :
PSNR value (dB) of the 8 McM images in Figure 4(a).

Table 3 :
PSNR value (dB) of the 8 Kodak images in Figure 4(b).