Joint Demosaicing and Denoising Based on Interchannel Nonlocal Mean Weighted Moving Least Squares Method

Nowadays, the sizes of pixel sensors in digital cameras are decreasing as the resolution of the image sensor increases. Due to the decreased size, the pixel sensors receive less light energy, which makes it more sensitive to thermal noise. Even a small amount of noise in the color filter array (CFA) can have a significant effect on the reconstruction of the color image, as two-thirds of the missing data would have to be reconstructed from noisy data; because of this, direct denoising would need to be performed on the raw CFA to obtain a high-resolution color image. In this paper, we propose an interchannel nonlocal weighted moving least square method for the noise removal of the raw CFA. The proposed method is our first attempt of applying a two dimensional (2-D) polynomial approximation to denoising the CFA. Previous works make use of 2-D linear or directional 1-D polynomial approximations. The reason that 2-D polynomial approximation methods have not been applied to this problem is the difficulty of the weight control in the 2-D polynomial approximation method, as a small amount of noise can have a large effect on the approximated 2-D shape. This makes CFA denoising more important, as the approximated 2-D shape has to be reconstructed from only one-third of the original data. To address this problem, we propose a method that reconstructs the approximated 2-D shapes corresponding to the RGB color channels based on the measure of the similarities of the patches directly on the CFA. By doing so, the interchannel information is incorporated into the denoising scheme, which results in a well-controlled and higher order of polynomial approximation of the color channels. Compared to other nonlocal-mean-based denoising methods, the proposed method uses an extra reproducing constraint, which guarantees a certain degree of the approximation order; therefore, the proposed method can reduce the number of false reconstruction artifacts that often occur in nonlocal-mean-based denoising methods. Experimental results demonstrate the performance of the proposed algorithm.


Introduction
Demosaicing refers to the task of reconstructing a full-color image from incomplete color samples sensed by the image sensor of a digital camera. The color filter array (CFA), which is placed over the image sensor, determines what colored pixels are missing. The most commonly used color filter array is the Bayer CFA [1]. One significant problem for the demosaicing process is the problem of noise. Nowadays, the sizes of the pixel sensors are decreasing with the increase in the resolutions of camera sensors. This makes the pixel sensors more sensitive to noise as the smaller sensors often receive less approximation for denoising of medical images. In [26], a modified singular value thresholding with minimizing error constraints is used for the segmentation of noisy signals by an orthogonal polynomial approximation.
In this paper, we propose a nonlocal interchannel weighted moving least square method for the denoising of noisy Bayer CFAs. The proposed method is a first attempt of applying a two-dimensional polynomial approximation-based denoising on noisy CFA patterns. Previous conventional joint denoising and demosaicing methods take either a 2-D approach with nonpolynomial approximation, e.g., a two-dimensional linear weighted sum of the neighboring pixels, or a directional 1-D polynomial approximation approach. The reason that previous polynomial approximation methods do not use a 2-D approach is due to the difficulty of the weight control in forming the 2-D approximation shape. If the conventional weighted least squares approach is applied to the 2-D polynomial approximation, the resulting 2-D approximation shape cannot represent the discontinuities of edge regions well, and therefore, will not result in a smooth image. To overcome this problem in our method, we compute the nonlocal mean weights for the weighted least square directly from the raw color filter array by taking the interchannel information into account. The goal of solving the weighted least squares problem is to find the filter coefficients for the 2-D filters that are used in the reconstruction of the pixel values; the better weights will result in finding better 2-D filters. After the filter coefficients are obtained, we again perform a nonlinear function on these coefficients so that the filtering of the data points with the 2-D filters only involves the important data points.
Compared with nonlocal-mean-based methods, e.g., the BM3D denoising, which only uses the mean filtering based on the nonlocal similarities of the pixels, the proposed method incorporates an extra reproducing constraint into the denoising scheme. Normally, nonlocal-mean-based denoising methods work well when there are iterative patterns in the image; however, if the correlations between the patches in the search region weaken, nonlocal-mean-based denoising methods can produce many false intensity values. This is because the denoising is performed without considering the approximation error. In comparison, the reproducing constraint is directly related to the reconstruction error, i.e., the approximation error. Using a reproducing constraint is better than nonlocal-mean-based methods in that it prevents the approximated denoised CFA image from being too different from the original CFA image. This improves the approximation accuracy.
We summarize the main contributions of the proposed method as follows: • For the first time, we applied a two-dimensional polynomial-approximation-based denoising on noisy CFA patterns.

•
Compared with nonlocal-mean-based methods, e,g., the BM3D denoising, the proposed method incorporates an extra reproducing constraint into the denoising scheme. This guarantees an approximation accuracy to a desired order.

•
We incorporate interchannel information into the polynomial approximation by determining the nonlocal weights directly from the noisy raw CFA image.

Relation of the Proposed Work to Sensors
Nowadays, most digital cameras acquire images with a single monochrome image sensor overlaid by a color filter array (CFA) to capture color information ( Figure 1). The reason for the use of a single image sensor in digital cameras is to reduce the cost. However, due to the use of a single image sensor, the color channels are undersampled and the missing color information has to be restored. The aim of a demosaicing algorithm in the IP (Image Processing) module in digital cameras is to reconstruct the full color image from the spatially undersampled color channels output from the CFA. As the original CFA image is noisy and of low resolution, demosaicing algorithms are centered on resolution improvement and denoising. This paper proposes a joint denoising and demosaicing framework based on the moving least square (MLS) method and, therefore, helps to overcome the physical limitations of single monochrome image sensors in digital cameras.

Residual Interpolation
The residual interpolation (RI) method proposed in [4] is a demosaicing method with excellent performance. It is an algorithm developed by integrating the residual interpolation method with the gradient-based threshold-free (GBTF) algorithm [6]. With this method, first, a tentative G pixel value is estimated at the positions of the R and B pixels by the Hamilton and Adams interpolation as in the GBTF algorithm. After that, the R and B pixel values are interpolated by the residual interpolation process. For example, for the R image, the RI algorithm first generates the tentative estimate (R) of the R image by a guided upsampling process, and then estimates the residuals between the R pixel and thẽ R pixel values at the R pixel positions. The residuals are also interpolated to result in the interpolated residual image. After that, the demosaiced R imageR is constructed by addingR to the interpolated residual image. The demosaiced B imageB is constructed in the same manner. In the proposed method, we utilize the RI method to construct the full-color image after the denoised CFA image has been obtained with the proposed interchannel nonlocal mean weighted moving least squares method.

Moving Least Square Methods with Total Variation Minimization
Moving least squares refers to the method of reconstructing a continuous function from a set of unorganized point samples by calculating a weighted least squares measure around the point at which the reconstruction is required. The solution of the MLS method has a closed-form and is easily computed by solving a linear system. It was shown to be quite useful in interpolation-based image-processing such as super-resolution and image zooming [27][28][29]. The MLS method has been applied to 2-D linear and nonlinear systems [30] for the interpolation in meshless environments [31] and for other image processing tasks such as nonlinear color transfer [32], but not to the problem of denoising the CFA image, since one of the major drawbacks of the MLS method is that it is weak against the noise. To overcome the problem of weakness against noise, in [33], we proposed the incorporation of the total variation regularization [34] into the MLS framework for better denoising power; however, if the MLS formulation in [33] is directly applied to the joint demosaicing and denoising problem, it yields poor results, as the correlation between the color channels cannot be well measured because two-thirds of the data are missing and the positions of the missing data are different for each color channel. Therefore, to apply the MLS method to denoise the color filter array, the formulation has to be changed to fit the problem and the weights have to be calculated while taking the interchannel information into account.

Problem Formulation
Let the ground-truth color image on an image domain Ω be given as The problem of demosaicing is to construct a color image u * from a noisy CFA, x, where x is the addition of the mosaiced pattern image u and the noise n: It should be noted that u is a color image, which is why we use a vector representation for u, while the mosaiced pattern image u is a monotone CFA image. In this paper, the mosaiced pattern image u is a Bayer pattern image and n is a white Gaussian noise which follows a normal distribution n(i, j) ∼ N (0, σ 2 ) with a certain variance σ. Our goal is to reconstruct u * as close as possible to the true image u. For the denoising of the CFA, we first construct a local polynomial approximation function for each pixel that well reflects the local structure at that pixel. We then take the coefficients of the local polynomial approximation function and apply a nonlinear transform on these coefficients. After that, we take the dot product of the filter with a neighborhood of pixels to decide the value of the centered pixel. The filtered pixels become the denoised version of the original noisy pixels in the CFA. We describe the proposed method in detail in the following subsection.

Interchannel Data Weighted Least Squares Reconstruction
Let x be a given mosaiced CFA image and let p C = (p C , q C ) denote a pixel position corresponding to the color C ∈ {R, G, B}. We want to construct the local polynomial approximation functions L p C (r) corresponding to each color C ∈ {R, G, B}: where Ω p C denotes the set of pixels that are in the neighborhood of p C and have the same color as p C . For example, Ω p G := Ω p C=G contains only the Green pixels in the CFA that are in the neighborhood of p C=G , i.e., Ω p G contains the nearest 41 green pixel positions (p, , q) which are (p, Then, the polynomial L p C is constructed to match the distribution of the data in the set Ω p C . By doing so, L p C takes into account the local structure at p C . For the Red and Blue pixels in the CFA, we construct the polynomials L p R and L p B from the data of larger regions, i.e., from the nearest 49 Red pixels or Blue pixels, respectively. For a specific Red pixel p R = (p R , q R ), the nearest 49 Red pixels are The polynomial L p G (p G ) can be obtained by the following minimization, which takes the total variation regularization into account: It should be noticed that the gradients ∇L p G (p G + 2i, q G + 2j) and ∇L p G (p G + 2i − 1, q G + 2j − 1) are computed by taking only the Green pixels in the CFA into account and not the neighboring pixels, as in the case of computing the gradients in normal images. For the computation of L p R (p R ) and L p B (p B ), we solve the following minimization problems, and argmin where once again the gradients ∇L p R (p R + 2i, q R + 2j) and ∇L p B (p B + 2i, q B + 2j) are computed with respect to the Red pixels and Blue pixels in the CFA, respectively. The weighting functions for Equation (1), and for Equations (2) and (3), respectively, where G σ (·, ·) denotes a Gaussian function with standard deviation σ. The nonlocal weight θ p C measures the similarity of the data structure of p C = (p C , q C ) andp C = (p C ,q C ), where p C = (p C , q C ) andp C = (p C ,q C ) are the pixels that belong to the same color; however, even though θ p C measures the similarity between pixels of the same color, all the RGB pixels are involved in the computation of the similarity. For example, even thoughp G = (p G ,q G ) is the position of the pixel belonging to the Green color, (i, j), i =p G − 3, . . . ,p G + 3, j =q G − 3, . . . ,q G + 3 are the pixel positions of all colors. This is in contrast with (1)- (3), where we used only the pixels corresponding to a specific color. The incorporation of all the RGB pixels provides the computation of the similarity with interchannel information. This enhances the accuracy in computing the weighting functions. As a different weighting function results in a different local polynomial approximation function, interchannel information results in a better polynomial approximation. The size of the window that contains the pixels involved in the computation is 7 × 7 for the computation of θ p G and 9 × 9 for the computation of θ p R and θ p B , as shown in Figure 2. Normally, after obtaining all the local polynomial approximation functions L p C (r) corresponding to all colors C ∈ {R, G, B} and all pixels r, we already can obtain the denoised CFA image. This can be done by evaluating the values of all the functions L p C at the points where the functions have been constructed, i.e., evaluate the values L p C (p C ) for all p C . By integrating all the values L p C (p C ) into a 2-D image, we obtain the denoised CFA image; however, with the proposed method, we take one more step to eliminate the effect of the data points that contribute little to the reconstruction of the function values L p C (p C ). This can be done by first utilizing the theorem that the L p C (p C ) can be expressed as a dot product between the filter weights and the data points [35]. Applying the theorem to our case, we can rewrite the function values By inserting (4)- (6) into (1) where The above transform filters out the filter coefficients that are smaller than a predefined threshold value th, and renormalize the remaining filter coefficients. By filtering out the small filter coefficients, the data points that have little contribution to the construction of L p C (p C ) are excluded from the construction process. This prevents unrelated pixels from affecting the construction process, which prevents an oversmoothed CFA. Measuring the local structure similarity using (a) 9 × 9 windows for R and B pixels and (b) 7 × 7 windows for G pixels. Now, using the transformed filter coefficients {φ even th (w G 2i,2j )}, {φ odd th (w G 2i−1,2j−1 )}, {φ th (w R i,j )}, and {φ th (w B i,j )} instead of the original weights, we construct L p C (p C ), C=R,G,B instead of L p C (p C ), C=R,G,B : Finally, we replace all the intensity values in {u(p C , q C )} p C ,q C ∈Ω, C=R,G,B in the noisy CFA with { L p C (p C , q C )} p C ,q C ∈Ω, C=R,G,B , where Ω is the domain corresponding to the CFA image. Here, (p C , q C ) refers to the pixel position corresponding to the color C. With the denoised CFA { L p C (p C , q C )} p C ,q C ∈Ω, C=R,G,B , we are now ready to apply the residual interpolation algorithm [4] to obtain the full-color image, i.e., the three-channel array u * = {u * (i, j)} i,j∈Ω , where u * (i, j) = [R * (i, j), G * (i, j), B * (i, j)] T is the full-color pixel reconstructed by the demosaicing method. Figure 3 shows the block diagram of the proposed joint denoising and demosaicing method, while Figure 4 shows how the input image is processed by the proposed method and also shows the intermediate image results at each stage.

Experimental Results
We compared the joint denoising and demosaicing results of the proposed scheme with the RI [4] to show the denoising power of the proposed scheme-the block-matching 3D (BM3D) filtering [7], which uses an enhanced sparse representation in transform-domain; and the alternating direction minimization multipliers (ADMM) method [9], which efficiently incorporates the total variation minimization and the BM3D filtering into a unified framework. The BM3D and the ADMM methods are state-of-the-art methods for denoising the color filter array. Both the ADMM and the BM3D use a denoising method that groups similar patches into a 3D volume, then performs a transform-domain shrinkage on this 3D volume. The collaborative filtering used in these methods can reveal even the finest details shared by grouped blocks, and therefore, can preserve essential and unique features of each block.
The experiments show that the proposed method is comparable to the state-of-the-art denoising methods and can better preserve some features that the ADMM and BM3D methods cannot. We experimented on the McMaster dataset, which contains images that are closer to the images taken with a real digital camera than the Kodak dataset. We used two noisy datasets of different noise levels, i.e., on a dataset with weak noise where the noise is generated from a zero-mean Gaussian distribution with σ = 7.65, and a dataset with a higher level of noise generated with σ = 12.75, where σ is the standard deviation.
The number of data points used in the reconstructions of the local polynomial approximation functions in the Green channel is 41, where the data are extracted from a 9 × 9 window. For the Red and the Blue channels, the number of data points is 49 for each channel, which is extracted from a 15 × 15 region. The degree of polynomials and the width of the Gaussian for the Green channel are 3 and 0.3, respectively, while for the Red and Blue channels, the degree of polynomials and the width of the Gaussian are both 3 and 0.5, respectively. Figure 5 shows a comparison of the denoising results on the MCM 9 image with weak noise (σ = 7.65). We see that even though the BM3D method has the highest PSNR value, there are some blurry artifacts, which can be observed in the enlarged images of Figure 5 and especially in the region corresponding to the fruit basket frame, as can be seen in Figure 5i. The reason that the BM3D results in blurry artifacts is due to the fact that the fruit basket image has little repeatable patterns, therefore, the 3D filtering in the BM3D blurs the thin frame regions. This is also true for the ADMM method, as can observed in Figure 5h. However, with the proposed method, the fruit basket frame is better reconstructed and more distinctive, as can be seen in Figure 5j, which is due to the fact that the proposed method applies the reproducing constraint in the reconstruction process. Tables 1-3 compare the PSNR, FSIM, and SSIM values between the different demosaicing methods. The tables show the PSNR, FSIM, and SSIM values for all 18 images in the McMaster dataset for different noise levels and different demosaicing methods. The proposed method is comparable to the state-of-the-art joint demosaicing and denoising methods. The BM3D method in [7] shows the largest PSNR, FSIM, or SSIM values for most images due to the large denoising power, but also shows some local artifacts when there is no repeatable pattern, or when the image structure is small. In this case, the proposed method shows more desirable results. This can be observed again in the dataset with more noise (σ = 12.75). Figures 6 and 7 show the results of the MCM 10 and 18 images. Again, it can be observed in the enlarged images of the ADMM (Figure 6h) and the BM3D (Figure 6i) methods that the small details become blurry, which can be observed especially inside the blue boxes. However, the small details are sharply reconstructed with the proposed method, as can be seen in Figure 6j. Figure 7 shows that the colors of the small-scaled structures are better preserved with the proposed method. It can be observed in the enlarged images of the ADMM (Figure 7h) and the BM3D (Figure 7i) methods that the brown color has been diffused to green, especially inside the blue boxes, whereas in the original image (Figure 7f) there are several small-scaled structures with brown colors. The small-scaled structures with brown colors are well preserved with the proposed method as can be seen in Figure 7j.

Conclusions
In this paper, we proposed a two-dimensional polynomial-approximation-based denoising method with nonlocal weights for the denoising of noisy color filter arrays (CFAs). It is difficult to construct a two-dimensional approximation function from the sparse data points in the CFA image that can represent the original color image well, since the missing information makes it difficult to extract the data points for the two-dimensional approximation that can reconstruct the full color image with high accuracy. The major contribution of this paper is the proposal of a method that can compute the weights that indicate the importance of data points with high accuracy, so that the data points that are close to the point at which the reconstructed value is requested can be extracted with high accuracy. To obtain a high-accuracy estimation of the data points and a high-accuracy two-dimensional polynomial approximation, we incorporated the interchannel information into the calculation of the similarity weights. This results in the reconstruction of digital color images with high visual quality. After the weights have been obtained, we applied a nonlinear transform on the filter coefficients and reconstructed the denoised CFA pixels. We also eliminated the data points that have little contribution to the construction of the local filters to exclude the unrelated pixels from the construction process, which prevents an oversmoothed CFA image. Due to the use of the interchannel nonlocal mean weights and the incorporation of the reproducing constraint and the transformed filtered coefficients, we could get an approximation function that effectively preserves the small-scaled features in the image that other conventional denoising schemes cannot preserve well. This is due to the fact that conventional denoising schemes do not possess the extra reproducing constraint, whereas with the proposed scheme, we could get a reproducing constraint with high order due to the highly accurate weights and data points. The high order of the proposed approximation method improves the resolution in the demosaiced image.
One of the advantages of the proposed method is that it can be easily extended to be used with other CFAs due to the characteristics of the MLS, which reconstructs a function from a set of unorganized point samples. Therefore, the proposed method can be well-combined with existing super-resolution or other interpolation methods. The use of the proposed method with other CFA formats will be one of our further research topics. Furthermore, if we use basis functions other than polynomial basis in the reproducing constraint term, there is much space left for improvement in the approximation accuracy.