4-Band Multispectral Images Demosaicking Combining LMMSE and Adaptive Kernel Regression Methods

In recent years, multispectral imaging systems are considerably expanding with a variety of multispectral demosaicking algorithms. The most crucial task is setting up an optimal multispectral demosaicking algorithm in order to reconstruct the image with less error from the raw image of a single sensor. In this paper, we presented a four-band multispectral filter array (MSFA) with the dominant blue band and a multispectral demosaicking algorithm that combines the linear minimum mean square error (LMMSE) and the adaptive kernel regression methods. To estimate the missing blue bands, we used the LMMSE algorithm and for the other spectral bands, the directional gradient method, which relies on the estimated blue bands. The adaptive kernel regression is then applied to each spectral band for their update without persistent artifacts. The experiment results demonstrate that our proposed method outperforms other existing approaches both visually and quantitatively in terms of peak signal-to-noise-ratio (PSNR), structural similarity index (SSIM) and root mean square error (RMSE).


Introduction
Digital color cameras generally sensitive to three bands of the visible electromagnetic spectrum are used to capture digital color images representing the reflectance of the observed object. Nowadays, technological advancement has made it possible to overcome this three-band limitation with the development of multispectral digital cameras to acquire multispectral images with more than three spectral bands per pixel. There are several types of multispectral image acquisition systems including single-sensor one-shot cameras which are equipped with a multispectral filter mosaic. However, in the raw image from the sensor, each pixel is characterized by a single available spectral band. We will have to reconstruct the missing spectral bands by the demosaicking method. The reconstruction performance depends on the optimal choice of MSFA and the multispectral demosaicking algorithm.
Several MSFA patterns are proposed in the literature. To our knowledge, Miao et al.
[1] are the first to propose a generic MSFA model from a binary tree by recursively separating the checkerboard pattern based on a tree decomposition which defines the number of spectral bands and the probability of occurrence of each band. Aggarwal et al. [2] meanwhile implemented two MSFA patterns, one random and the other uniform, which can be generalized to any number of bands. In [3], Monno et al. proposed a five-band MSFA based on the dominant G band requirement, which is used by Jaiswal et al. in their multispectral demosaicking algorithm [4]. To overcome the difficulties in combining spectral resolution and spatial correlation, Mihoubi et al. proposed a 16-band MSFA without a dominant spectral band [5]. Recently, Bangyong et al. designed a uniform four-band MSFA pattern [6] with the same probability of occurrence for each band and a nine-band MSFA pattern [7] in which one band is dominant and the other eight have the same probability of occurrence arranged in a 4 × 4 mosaic.
Many multispectral demosaicking algorithms using the designed MSFAs have been implemented in the literature [8]. Miao et al. [9] proposed a binary tree-based edgesensing (BTES) multispectral demosaicking algorithm that recursively performs binary tree-based edge detection interpolation. However, the performance of this algorithm in classical edge detection interpolation is limited. Recently, Monno et al. [10][11][12] proposed a series of demosaicking algorithms for its proposed five-band MSFA. The first of these algorithms [10] developed several guide images that were used in the interpolation of the different spectral bands. The authors used residual interpolation to generate a guide image for structure-preserving interpolation [11] and proposed adaptive residual interpolation by adaptively combining two algorithms based on residual interpolation and selecting an appropriate number of iterations for each pixel [12]. Jaiswal et al. [4] used the high-frequency component of the G-band to interpolate the other bands based on an inter-band correlation analysis while Mihoubi et al. [2] proposed a 16-band MSFA algorithm based on a pseudo-panchromatic image (PPI), which is estimated by applying an averaging filter to the raw image and then adjusted such that the PPI values are correlated. The difference between each available value of the adjusted raw image and PPI is calculated. The calculated local directional weights are then used to estimate the fully defined difference using an adaptive weighted bilinear interpolation. Each band is finally estimated by adding a PPI and the difference. In [6], a method of applying directional interpolation along the edges of an image was proposed. In this method, the image edges are calculated from the raw image to define the direction interpolation with the neighbors. Considering the features of the filter arrays, image edges, and a constant hue, the missing bands per pixel were recovered from the existing bands. Then, the image is separated into high-and low-frequency components by applying a wavelet transform, and the high-frequency images that are highly correlated are modified using luminance information to refine the demosaicked image. In [7], a multispectral algorithm that estimates the missing dominant band at each spatial position with a weighted average of the neighboring values of the dominant band was described. The dominant band reconstructed at different spatial positions is then used as a guided image to estimate all other missing bands using the guided filter and residual interpolation.
Multispectral images demosaicked using the previous algorithms suffer from severe artefacts in edge regions. To overcome these limitations, a new avenue of multispectral demosaicking called the LMMSE method is being explored. Zhang and Wu [13], in the demosaicking of color images such as Bayer's CFA [14], had developed the LMMSE method which is based on the assumption that the gradient of the G and R/B channels correspond to low-pass filtering, given their strong correlation. The LMMSE adaptively estimates the missing G values in both horizontal and vertical directions and then merges them optimally. A very interesting result is the introduction of the neighborhood in the LMMSE formulation by Amba et al. [15] for color images they recently extended to eight-band multispectral demosaicking by applying a linear operator that minimizes the root mean square error between the reconstructed image and the original raw image [16]. This linear operator multiplied by the MSFA image provides an estimate of the reconstructed image. According to [16], the LMMSE method constitutes a good potential candidate for real-time applications because, after training, it could be integrated into the equipment of the camera and operate in real time without losing the generality required by the various provisions present on the market.
The contributions of our paper focused on the LMMSE method in [13] and the adaptive kernel regression kernel as described in [17] are of three kinds. (1) We identified with justification a generic four-band MSFA with the dominant blue band for our multispectral demosaicking algorithm. (2) We proposed the directional LMMSE method for estimating the missing blue bands and the directional gradient method for the other three spectral bands. (3) To take into account the details at the edges and the denoising of the reconstructed image, we have successfully combined the LMMSE method with the adaptive kernel regression. This paper is organized as follows: in the second section, we justify our proposed four-band MSFA and the application assumptions of the LMMSE method. The existing LMMSE method and the adaptive kernel regression used are described in the third section. The proposed algorithm and the experimental results are, respectively, presented in the fourth and fifth sections. based on the binary tree with the probability of occurrence of each spectral band ( Figure 1). The multispectral images of the cave dataset [18] used in our simulations are acquired with a camera whose sensor is fitted with the liquid crystal tunable filter (LCTF) [19] (Figure 2) such as the energy of the blue band of wavelength λ = 450 nm is very weak compared to the other bands, followed by the orange band whose wavelength is λ = 600 nm. The red band (λ = 700 nm) has the greatest energy preceded by the green band (λ = 550 nm). According to [20], the energy imbalance between the different spectral bands produces in the demosaicked image, severe degradation of the lowenergy bands due to their sensitivity to noise. It then appears necessary that a balancing be carried out to optimize the shape of the transmittance filters. To balance the energies of the different spectral bands and avoid degradation of the low energy bands, we opted for an MSFA with a dominant blue band with a probability of occurrence of 1/2 followed by the orange band with a probability of occurrence of 1/4, and 1/8 is the probability of occurrence of the red and green bands. reconstructed image, we have successfully combined the LMMSE method with the adaptive kernel regression. This paper is organized as follows: in the second section, we justify our proposed four-band MSFA and the application assumptions of the LMMSE method. The existing LMMSE method and the adaptive kernel regression used are described in the third section. The proposed algorithm and the experimental results are, respectively, presented in the fourth and fifth sections.

Design of the 4-Band MSFA
The identified four-band MSFA is based on the generic method of Miao et al.
[1] based on the binary tree with the probability of occurrence of each spectral band ( Figure  1). The multispectral images of the cave dataset [18] used in our simulations are acquired with a camera whose sensor is fitted with the liquid crystal tunable filter (LCTF) [19] (Figure 2) such as the energy of the blue band of wavelength λ = 450 nm is very weak compared to the other bands, followed by the orange band whose wavelength is λ = 600 nm. The red band (λ = 700 nm) has the greatest energy preceded by the green band (λ = 550 nm). According to [20], the energy imbalance between the different spectral bands produces in the demosaicked image, severe degradation of the low-energy bands due to their sensitivity to noise. It then appears necessary that a balancing be carried out to optimize the shape of the transmittance filters. To balance the energies of the different spectral bands and avoid degradation of the low energy bands, we opted for an MSFA with a dominant blue band with a probability of occurrence of 1/2 followed by the orange band with a probability of occurrence of 1/4, and 1/8 is the probability of occurrence of the red and green bands. After balancing the energies of the different spectral bands according to the proposed MSFA pattern, their spectral sensitivities are shown in Figure 3. After balancing the energies of the different spectral bands according to the proposed MSFA pattern, their spectral sensitivities are shown in Figure 3.

Application Assumptions of LMMSE Method
The LMMSE method as used for RGB images obeys two assumptions. Firstly, in natural images, the different spectral bands are strongly correlated. Then, the gradient of the different bands remains constant and constitutes a smooth process (low pass) [13]. To verify these assumptions, we first determined in Table 1, the spectral correlation of the different spectral bands of the cave dataset [18] multispectral images used. The spectral correlation of two bands is best if the correlation coefficient between these two bands is between 0.5 and 1. From the analysis of Table 1 and with a few minimal exceptions, all the spectral bands of the different multispectral images are strongly correlated. Secondly, we have plotted the power spectral functions of the gradients of the different bands for three multispectral images (Figures 4-6). When we analyze these different spectral functions, we have realized that the power of the different gradient signals is concentrated in the low-frequency band then that each of these functions has a peak around the zero frequency. Thus, the two above-mentioned are therefore verified for our simulation multispectral images.

Application Assumptions of LMMSE Method
The LMMSE method as used for RGB images obeys two assumptions. Firstly, in natural images, the different spectral bands are strongly correlated. Then, the gradient of the different bands remains constant and constitutes a smooth process (low pass) [13]. To verify these assumptions, we first determined in Table 1, the spectral correlation of the different spectral bands of the cave dataset [18] multispectral images used. The spectral correlation of two bands is best if the correlation coefficient between these two bands is between 0.5 and 1. From the analysis of Table 1 and with a few minimal exceptions, all the spectral bands of the different multispectral images are strongly correlated. Secondly, we have plotted the power spectral functions of the gradients of the different bands for three multispectral images (Figures 4-6). When we analyze these different spectral functions, we have realized that the power of the different gradient signals is concentrated in the low-frequency band then that each of these functions has a peak around the zero frequency. Thus, the two above-mentioned are therefore verified for our simulation multispectral images.

Application Assumptions of LMMSE Method
The LMMSE method as used for RGB images obeys two assu ural images, the different spectral bands are strongly correlated. T different bands remains constant and constitutes a smooth process ify these assumptions, we first determined in Table 1, the spectra ferent spectral bands of the cave dataset [18] multispectral images  Table 1. Spectral correlation coefficients for red/green (C rg ), red/blue (C rb ), red/orange (C ro ), blue/green (C bg ), blue/orange (C bo ) and orange/green (C og ).

LMMSE Demosaicking Method
In their color image demosaicking algorithm, Zhang and Wu [13] used the LMMSE method to estimate the missing G bands at different pixels. We briefly present here the estimate of the G-band at each spatial position , of the red pixels in the CFA image. Thus, the missing G-band at the red pixels is obtained according to the formula used: , , The gradient ∆ , of the red and green bands is estimated in both horizontal and vertical directions such as:

LMMSE Demosaicking Method
In their color image demosaicking algorithm, Zhang and Wu [13] used the LMMSE method to estimate the missing G bands at different pixels. We briefly present here the estimate of the G-band at each spatial position , of the red pixels in the CFA image. Thus, the missing G-band at the red pixels is obtained according to the formula used: , , The gradient ∆ , of the red and green bands is estimated in both horizontal and vertical directions such as:

LMMSE Demosaicking Method
In their color image demosaicking algorithm, Zhang and Wu [13] used the LMMSE method to estimate the missing G bands at different pixels. We briefly present here the estimate of the G-band at each spatial position (i, j) of the red pixels in the CFA image. Thus, the missing G-band at the red pixels is obtained according to the formula used: The gradient ∆ g,r of the red and green bands is estimated in both horizontal and vertical directions such as: Before calculating the gradient, a second order Laplacian interpolation is used beforehand to know at each pixel, the "missing samples". The noises associated with the directional gradient estimate are determined as: So, we have: The gradient ∆ g,r is estimated by the LMMSE method. Let us denote by x the gradient ∆ g,r , y the associates∆ g,r h and∆ g,r v and ϑ the associated noises ε g,r h and ε g,r v . The Equation (5) becomes: The optimal estimate of the minimum mean square error (MMSE) of x is defined: However, in practice, the probability p(x/y) is rarely known making estimation of MMSE difficult. Therefore, instead of MMSE, the authors used the LMMSE method to estimate x, such as:x E[x] is the mathematical expectation of x, Cov(x, y) the covariance and Var(y) the variance of y. By setting µ x = E[x], σ x 2 = Var(x) and σ ϑ 2 = Var(ϑ), the Equation (8) becomes: For an optimal estimate ofx, the latter is estimated adaptively by merging the values determined in both horizontal and vertical directions in the neighborhood of y. Denoting byx h (n) andx v (n) the horizontal and vertical LMMSE estimates of x obtained from the Equation (9), then by w h and w v the both horizontal and vertical weights, respectively, the optimal LMMSE estimate of x is defined by: With w h (i, j) + w v (i, j) = 1 to minimize the estimation error. x h and x v are the estimation errors ofx h andx v such that: σ x h 2 and σ x v 2 are, respectively, the variances of x h and x v . More information is given in [13].

Kernel Regression Method
Takeda et al. [17] proposed a kernel regression that is used in the iterative reconstruction of color images, and which takes into account limitations such as strong denoising along the edges, the high retention of detail in the edges, and the limited presence of blur in the reconstructed image. The estimate of y pixel at x i location is defined as: ∈ i is the associated noise and z(.) the regression function obtained by Taylor expansion of N-order.
vech(.) is a half-vectorization operator of the lower triangular portion of a symmetric matrix such as: The β n are obtained as below: They are computed by the following optimization problem: K is the kernel function and H the 2 × 2 smoothing matrix of order defined by: h is a global smoothing parameter, µ i a local density parameter which controls the kernel size and I an identity matrix.

Adaptive Kernel Regression
The adaptive kernel regression is an extension of the classical kernel regression [17] and structured in the same way as in Equation (17) where the classical kernel is replaced by the adaptive kernel.
H s = h s I is the spatial smoothing matrix. To avoid computational complexity, the order estimation is limited to N = 0. The necessary calculations are then limited to those which estimate the parameter β o such as: The value of a spectral band at a spatial position is determined by: Expressing K adapt in spatial and radiometric terms weakens the performance of the estimate. Consequently, the adaptive kernel is replaced by an adaptive steering kernel, the denoising of which takes place most strongly along the edges.
The steering matrix is defined as: The C i are symmetric covariance matrix used to temper the blurring effect around edges and whose values are obtained by a differentiation between the value of the central pixel and those of the neighboring pixels. The global smoothing parameter h makes it possible to have a strong denoising effect and the steering kernel is a Gaussian kernel.
where z x 1 (·) and z x 2 (·) are the first derivatives along x 1 and x 2 directions and w i is a local analysis window around the position of interest. We set the smoothing parameter h i to 2 to have a strong denoising effect along edges and the local density parameter µ i to 1 for kernel size control.

Proposed Multispectral Demosaicking Method
The proposed algorithm is subdivided into six main steps. The blue band being the dominant band of our MSFA (Figure 1c), this band is the first one estimated at the other pixels.

Blue Band Estimation by LMMSE Method
We estimate the blue band missing at the orange pixel by applying the formula: The gradient ∆ b,o is interpolated by the LMMSE method of Equations (10) and (11). We adopt the same strategy to estimate blue bands at red and green pixels such as: in G pixels (28)

Orange Band Estimation at Red and Green Pixels
The green and red bands have identical neighborhoods.
The gradient values ∆ b,o in the four directions ( Figure 7) northwest (nw), northeast (ne), southwest (sw) and southeast (se) in the neighborhood of a green or red pixel are, respectively, denoted by:

Orange Band Estimation at Red and Green Pixels
The green and red bands have identical neighborhoods.
The gradient values ∆ , in the four directions (Figure 7) northwest , northeast , southwest (sw) and southeast in the neighborhood of a green or red pixel are, respectively, denoted by: Each orange band is then estimated at the green and red pixels by the formula:

Green Band Estimation at Red Pixels and Vice Versa
We apply the same strategies as before but in a wider neighborhood (Figure 8) in the north (n), south (s), east (e), and west (w) directions.  Each orange band is then estimated at the green and red pixels by the formula:

Green Band Estimation at Red Pixels and Vice Versa
We apply the same strategies as before but in a wider neighborhood (Figure 8) in the north (n), south (s), east (e), and west (w) directions.
The gradient ∆ , is interpolated by the LMMSE method of Equations (10) and (11). We adopt the same strategy to estimate blue bands at red and green pixels such as:

Orange Band Estimation at Red and Green Pixels
The green and red bands have identical neighborhoods.
The gradient values ∆ , in the four directions (Figure 7) northwest , northeast , southwest (sw) and southeast in the neighborhood of a green or red pixel are, respectively, denoted by: Each orange band is then estimated at the green and red pixels by the formula:

Green Band Estimation at Red Pixels and Vice Versa
We apply the same strategies as before but in a wider neighborhood (Figure 8) in the north (n), south (s), east (e), and west (w) directions.

Red and Green Bands Estimation at Orange Pixels
The Figure 9 shows the neighborhood of orange band for original and estimated red and green pixels.

Red and Green Bands Estimation at Orange Pixels
The Figure 9 shows the neighborhood of orange band for original and estimated red and green pixels. Taking into account the neighborhood of the orange bands for the red and green pixels (Figure 9), original or estimated, we estimate the red and green spectral bands for the orange pixels according to the formulas: , , ∆ ,

Red, Green, and Orange Bands Estimation at Blue Pixels
From Figure 10b-d, we can see a symmetry of the neighborhood of the blue pixels for the red, green and orange pixels. If we denote by ∆ ; ∆ ; ∆ and ∆ the directional gradients of the blue and orange pixels in a neighborhood, we compute the average of the gradient bilinearly as follows: Therefore, the missing orange band at the blue pixels is estimated by the relation: Similarly, we estimate the red and green bands at the blue pixels such as: , , ∆ , Taking into account the neighborhood of the orange bands for the red and green pixels (Figure 9), original or estimated, we estimate the red and green spectral bands for the orange pixels according to the formulas: where∆ ij and∆ bg are respectively the gradient values in the four directions (Figure 9b,c) for the red and green pixels.

Red, Green, and Orange Bands Estimation at Blue Pixels
From Figure 10b-d, we can see a symmetry of the neighborhood of the blue pixels for the red, green and orange pixels. If we denote by ∆ n bo ; ∆ s bo ; ∆ w bo and ∆ e bo the directional gradients of the blue and orange pixels in a neighborhood, we compute the average of the gradient bilinearly as follows:

Estimated Bands Enhancement Using Adaptive Kernel Regression
Although the previous estimation formulas have worked well for color images as in [13], their use for multispectral images remains limited, especially not taking into account details in the strong edge or rich texture. To correct these imperfections, each estimated spectral band is refined by using the adaptive kernel where we replaced the refined Equation (21) with the formula as defined in [21]. Thus, the refined a spectral band at location is defined by: where is the set of neighbor pixel locations of the location , the sampled value at the location , the binary mask at the location that set to one if the data are sampled at an associated location and set to zero otherwise and is the nor- Therefore, the missing orange band at the blue pixels is estimated by the relation: Similarly, we estimate the red and green bands at the blue pixels such as:

Estimated Bands Enhancement Using Adaptive Kernel Regression
Although the previous estimation formulas have worked well for color images as in [13], their use for multispectral images remains limited, especially not taking into account details in the strong edge or rich texture. To correct these imperfections, each estimated spectral band is refined by using the adaptive kernel where we replaced the refined Equation (21) with the formula as defined in [21]. Thus, the refined a spectral band at x p location is defined by: where N p is the set of neighbor pixel locations of the location x p , S(x i ) the sampled value at the location x i , M(x i ) the binary mask at the location x i that set to one if the data are sampled at an associated location and set to zero otherwise and w x p is the normalizing factor, which is the sum of kernel weights. The adaptive steering kernel K H i steer is computed according to Equation (25) and the covariance matrix C x p according to Equation (26).

Experimental Results
In our experiments, we used 26 images from the 32 (the others being resemblances) of the cave dataset [18], in which multispectral images consist of 31-band multispectral images acquired under illuminant D65. The 31-band images were acquired every 10 nm at between 400 and 700 nm. The image size was 512 × 512 pixels. The CAVE dataset is often used as a standard multispectral image dataset.
To evaluate the performance of the proposed algorithm, we compared it with recent four-band multispectral demosaicking methods, namely generic binary tree edge sensing (BTES) [9], directional filtering and wavelet transformation (DFWF) [6], adaptive spectralcorrelation based demosaicking (ASCD) [4] and neighborhood in linear minimum mean square error (N-LMMSE) [16]. ASCD and N-LMMSE are a five-band and eight-band methods, respectively, which we implemented to four-band for comparison purposes. Visual and objective evaluations were also conducted.

Visual Performance Evaluations
For evaluation purposes, we selected four images with detailed structures as shown in Figures 11-14. From the partially zoomed-in view images (red areas in original images), one notes the visible presence of blurring and false colors artifacts in the images demosaicked with the algorithms BTES, DFWF, ASCD, and N-LMMSE as is the case in Figures 12b-e, 13b-e and 14b-e which, respectively, display the green, blue, and orange bands of the feathers, hairs, and cloth images. In Figure 11b-e showing the red band of the face image, these artifacts are more visible with the BETS and N-LMMSE algorithms. In Figure 12, we note the presence of ghost noise in part of the reconstructed images with the BTES, DFWF, ASCD, and N-LMMSE algorithms. The quality of the reconstructed image is considerably reduced by these artifacts which are due to the lack of edge-preserving of the BTES, DFWF, ASCD, and N-LMMSE algorithms. Our proposed method reconstructs images without significant blurring or zipper artifacts (Figures 11f, 12f, 13f and 14f). The four reconstructed images with our proposed demosaicking algorithm preserve details at edges and in textured areas better than the other algorithms. Overall, by comparing the results of the visual assessment, we can confidently say that our proposed method is better than the BTES, DFWF, ASCD, and N-LMMSE algorithms. reduced by these artifacts which are due to the lack of edge-preserving of the BTES, DFWF, ASCD, and N-LMMSE algorithms. Our proposed method reconstructs images without significant blurring or zipper artifacts (Figures 11f-14f). The four reconstructed images with our proposed demosaicking algorithm preserve details at edges and in textured areas better than the other algorithms. Overall, by comparing the results of the visual assessment, we can confidently say that our proposed method is better than the BTES, DFWF, ASCD, and N-LMMSE algorithms.   DFWF, ASCD, and N-LMMSE algorithms. Our proposed method reconstructs images without significant blurring or zipper artifacts (Figures 11f-14f). The four reconstructed images with our proposed demosaicking algorithm preserve details at edges and in textured areas better than the other algorithms. Overall, by comparing the results of the visual assessment, we can confidently say that our proposed method is better than the BTES, DFWF, ASCD, and N-LMMSE algorithms.

Quantitative Performance Evaluations
To quantitatively assess the objective performance of our proposed algorithm, we used the PSNR, SSIM, and RMSE metrics as described in [6,22] and calculated from the original and demosaicked images. The average values of PSNR, SSIM, and RMSE obtained from various algorithms are shown in Tables 2-4, respectively, such that the best scores are in bold. Note that the lower the value of RMSE, the better the performance of the algorithm.

Quantitative Performance Evaluations
To quantitatively assess the objective performance of our proposed algorithm, we used the PSNR, SSIM, and RMSE metrics as described in [6,22] and calculated from the original and demosaicked images. The average values of PSNR, SSIM, and RMSE obtained from various algorithms are shown in Tables 2-4, respectively, such that the best scores are in bold. Note that the lower the value of RMSE, the better the performance of the algorithm.

Quantitative Performance Evaluations
To quantitatively assess the objective performance of our proposed algorithm, we used the PSNR, SSIM, and RMSE metrics as described in [6,22] and calculated from the original and demosaicked images. The average values of PSNR, SSIM, and RMSE obtained from various algorithms are shown in Tables 2-4, respectively, such that the best scores are in bold. Note that the lower the value of RMSE, the better the performance of the algorithm. A careful analysis of Table 2 shows that the DFWF algorithm produced the highest PSNR value which is 46.9371 while it is 42.3519 for the proposed method. Moreover, this method gives better PSNR scores for twelve images out of twenty-six while ours produces ten images. The N-LMMSE method gives no good score but the other two algorithms for two images each. Images for which the PSNR values are high for the DFWF algorithm are smoother images, but not textured images. The results of the average values of the PSNR show that our proposed method ranks first in the competition with DFWF, while N-LMMSE comes in last. In Table 3, according to the SSIM values, our algorithm outperforms all others with a higher value of the metric for both sixteen out of twenty-six images and for the average value. It is followed by BTES with a better score of the SSIM for ten images, the other three methods come in last. In Table 4, our algorithm produced better RMSE scores for seventeen images, unlike the others which have low scores. The RMSE values confirm the previous results with a lower mean value for our proposed method than the other algorithms. We can therefore see that overall, the values of PSNR, SSIM, and RMSE obtained with our algorithm are better than those of the BTES, DFWF, ASCD, and N-LMMSE methods.

Conclusions
In this study, we identified a four-band MSFA pattern for single-sensor cameras, arranged in a 6 × 6 moxel half filled with the blue band taking into account the properties of the liquid crystal tunable filter with which the sensor surface of the camera used is covered to acquire simulation images of the cave dataset. Based on the existing one, we then proposed a consequent algorithm that combines the LMMSE method and the adaptive kernel regression. In the proposed algorithm, we estimated the missing blue bands by the LMMSE method and the other spectral bands by the directional gradient method which relies on the estimated blue bands. Finally, applying the adaptive kernel regression gradient method to each spectral band refines the band by ridding it of artifacts that can adversely affect the reconstruction performance. In the experiment, we evaluated the proposed algorithm both visually and quantitatively with the existing algorithms BTES, DFWF, ASCD, and N-LMMSE. The results show that our proposed algorithm outperforms the others both visually and in terms of PSNR, SSIM, and RMSE.
The future work consists in deepening the algorithm in terms of the number of spectral bands by varying the moxels such as they are of 4 × 4 and 8 × 8 support to fully appreciate the degradation behavior of each spectral band.