Directional Interpolation Based on Improved Adaptive Residual Interpolation for Image Demosaicking

1. Introduction

To save manufacturing costs and reduce the size of the device, most digital imaging devices use monochrome sensors that can only obtain the change of single color. Therefore, an algorithm must be used to estimate the loss of color to obtain a full-color image. This process is called demosaicking, which belongs to image interpolation. Image demosaicking is the key to imaging process, so it is especially necessary to improve the image quality by improving performance of the algorithm [1].

Since the image sensor captures the change of single color, Bayer [2] covered the photosensitive unit with an optical filter called color filter array (CFA). Currently, Bayer pattern CFA becomes the most commonly used model, and its schematic diagram is 5×5 pixels, which is shown in Fig. 1, where R, G, and B respectively represent color components obtained after passing through a filter of red, green, and blue. Their positions are marked as [TeX:] $$\mathrm{L}_{\mathrm{R}}, \mathrm{L}_{\mathrm{G}},$$ and [TeX:] $$\mathrm{L}_{\mathrm{B}},$$, respectively.

Based on Bayer pattern CFA algorithm, the earliest and most widely used interpolation algorithm was bilinear interpolation [3], which decomposes CFA into three mutually independent images and processes them independently. This algorithm is not complicated, which can suppress image overlapping effectively. However, in the area with complex edges, color artifacts and serrated effects are produced by bilinear interpolation [3], therefore, the quality of the image is usually imperfect. To improve this problem, researchers began to study the correlation of spectrum among different color components. Gunturk et al. [4] found that the correlation among different color components was mainly determined by smoothness of chromatic aberration, which means, the smoother the color plane, the better the image interpolation effect. This theory led to some algorithms based on color interpolation (CI) [5,6] being proposed. Zhang and Xu [5] proposed a directional linear minimum mean square-error (DLMMSE) estimation algorithm, which estimates the missing G component adaptively in different directions of the image. Although the DLMMSE algorithm obtains optimal color-difference evaluation at each individual pixel location, it ignores the additional information that may be provided by neighboring pixel locations. To obtain such additional information, Pkkucuksen and Altunbasak [6] proposed a gradient based threshold free (GBTF) algorithm, which uses a color-difference gradient instead of a color-difference variance to estimate the missing G component. The GBTF algorithm not only obtains additional information of neighboring pixel locations, but also eliminates the hard rules of DLMMSE, so the threshold setting is no longer needed. The GBTF algorithm assumes that the image has hyperspectral correlation. However, in natural environment, spectral correlation is usually weaker in the edge regions of the image than in the smooth regions.

In recent years, Kiku et al. [7] proposed the residual interpolation (RI), which uses a smoother residual plane than color plane for interpolation, and achieves a remarkable image restoration effect. Based on RI method, Kiku et al. [8] found that the smaller the Laplacian energy (LE) of the image, the better the image interpolation effect. Therefore, they proposed the minimized-Laplacian RI (MLRI) to minimize the value of LE. However, in these algorithms [7-8], the prediction error generated during the generation of G component affects the generation of R and B components. To improve this problem, Ye and Ma [9] proposed iterative RI (IRI) that tested the test datasets iteratively to reduce the estimation error. Kim and Jeong [10] combined the estimated pixel values with the inverse gradient weights in four directions, so the prediction error can be reduced. Yu et al. [11] weakened the effects of gradient weights in horizontal (H) and vertical (V) directions on the prediction error. Monno et al. [12] used the iterative thought of IRI to combine RI with MLRI, and then proposed adaptive RI (ARI) to generate G component, which improved the quality of image restoration. However, in the ARI algorithm, a non-adaptive MLRI algorithm is still adopted to obtain R and B components. Recently, Ni et al. [13] proposed a progressive collaborative representation (PCR) framework, which can incorporate any existing color image demosaicking method and further improve its demosaicking performance. Sun et al. [14] proposed a hybrid demosaicking algorithm based on fuzzy edge strength and RI method, which has better performance in terms of structural similarity and visual comparison.

To solve the problem of ARI algorithm, the adaptability of the algorithm is not only applicable to G component, but also to R and B components. In addition, the interpolation of R and B components includes not only H and V directions but also the diagonal direction. Meanwhile, it is more difficult to use the ARI algorithm directly for both types of interpolation directions simultaneously. According to disadvantages mentioned above, this paper proposes a solution. This paper takes the RI algorithm as research object. The main research contents include: understanding the principle of obtaining smooth residual plane by RI, MLRI, and IRI algorithms; reproducing ARI algorithm; improving the interpolation process of R and B components in the improved ARI (IARI) algorithm; testing and evaluating the image restoration effects of the five algorithms.

The rest of this paper is organized as follows. Section 2 starts with a brief review of RI algorithm basis and its derivative algorithm, and then ARI algorithm is described. The proposed method is designed in Section 3. Experimental results of the proposed method are presented in Section 4. Conclusions and remarks on possible further work are given finally in Section 5.

2.1 Residual Interpolation Algorithm Basis and Its Derivative Algorithm

The basic flow of RI algorithm is described in detail in [7]. The important details in RI algorithm include the interpolation of G component by the GBTF interpolation algorithm [6], and the guided filter (GF) of the image generated by G component after interpolation [15], where the GBTF interpolation algorithm is implemented on the basis of the Hamilton interpolation formula [15]. From [16], in fact, in the RI algorithm, the cost function minimizes the value of the residual. Therefore, the essence of the RI algorithm is to obtain the smooth residual plane by minimizing the residual. After that, based on the RI algorithm, Kiku et al. [8] proposed the MLRI algorithm according to the relation that the image smoothness index LE is negatively correlated with the RI effect. The basic flow of MLRI algorithm is described in detail in [8]. Moreover, it can be seen that in the MLRI algorithm, the cost function minimizes the image smoothness index LE, which is a measure of image smoothness. Therefore, the essence of the MLRI algorithm is to obtain the smooth residual plane by minimizing the image smoothness.

Later, Ye and Ma [9] proposed the IRI algorithm, which makes the residual plane smooth by iteration. The basic flow of IRI algorithm is described in detail in [9]. It can be seen that the residual plane becomes smoother and smoother during the iterative repetition. Therefore, the essence of the IRI algorithm is to minimize the residual by iterating RI’s minimized residual to obtain the smooth residual plane. Since the local window size gradually increases during iterative process, it is impossible to judge the iterative coefficient [17]. Hence, Ye and Ma [9-17] determined the local window size which is an important parameter of GF and has a huge effect on the performance of the algorithm, thereby determining the optimal iterative coefficient.

2.2 Adaptive Residual Interpolation Algorithm

The different principles of RI, MLRI, and IRI algorithms for obtaining smooth residual planes have advantages in images with different characteristics. Based on RI, MLRI, and IRI algorithms, the optimal algorithm at different pixel locations is adopted adaptively in ARI algorithm to obtain better G component interpolation images [12].

Since the calculation of G component at LR and LB is the same, this paper only introduces the calculation of G component at LR. The following part describes three steps in G component interpolation process of the ARI algorithm [12]. Since the interpolation and selection process of the iterative coefficients in H and V directions are the same, this paper only introduces the algorithm and selection process in the H direction.

Step 1: To obtain a linear G component, the RI and MLRI algorithms are utilized iteratively in H and V directions to obtain a linear G component. For each interpolated result, a set of directional interpolated G component images will be generated, and each of them corresponds to one iteration [12].

Step 2: For the interpolation results in each direction, the optimal iteration number is adaptively selected from the set of directional interpolated G component images at each pixel position, based on iterative rule proposed in [17]. The iterative criteria at the pixel location in the kth iteration can be expressed as [12]:

where the superscript [TeX:] $$(\cdot)^{\mathrm{H}}$$ represents the H direction, and the subscript represents the kth interpolation at [TeX:] $$(i, j)$$. Meanwhile, [TeX:] $$\boldsymbol{d}_{(i, j), k}^{\mathrm{H}}$$ represents the difference between the estimated pixel value and the previous interpolation results, and [TeX:] $$\delta d_{(i, j), k}^{\mathrm{H}}$$ represents the smoothness of the difference. They can be denoted as [12]:

where the superscripts [TeX:] $$(\cdot)^{\mathrm{GH}}$$ and [TeX:] $$(\cdot)^{\mathrm{RH}}$$ represent the residual calculation process of R and B components in H direction, respectively. In addition, the parameters [TeX:] $$(m, n)$$ are set to the optimal solution obtained in [17]. The optimal iteration coefficient [TeX:] $$k_{\text {best }}$$ for each pixel location can be adaptively expressed as [12]:

where the function [TeX:] $$g(\cdot)$$ represents spatial Gaussian smoothing filter with the standard deviation [TeX:] $$\sigma=2$$ [12]. When [TeX:] $$\boldsymbol{k}_{\text {best }}$$ is determined, according to the interpolation results of Step 1, the optimal results can be obtained by interpolating RI or MLRI algorithm corresponding to each pixel position in the H or V direction.

Step 3: Weighted average at each LR can adaptively combine the results of RI and MLRI algorithms interpolated in different directions, and the final interpolated image [TeX:] $$\tilde{\boldsymbol{G}}$$ will be obtained.

3. Improved Adaptive Residual Interpolation Algorithm

The adaptability of the conventional ARI algorithm is only applied to the interpolation process of G component, and the interpolation of R and B components still uses the conventional MLRI algorithm. When the conventional ARI algorithm performs iterative interpolation on G component, it only contains interpolation directions of H and V. However, the interpolation of R and B components includes not only H and V directions, but also the diagonal direction [12]. Moreover, it is difficult to use the ARI algorithm directly for two types of interpolation directions simultaneously. To solve this problem, the interpolation process of R and B components is divided into two categories according to the direction types [12]: (i) H and V directions and (ii) diagonal direction. The ARI algorithm is performed once in two categories, respectively. The processing of B component is the same as R component, and then the interpolation process of R component includes the following two steps:

Step 1: The ARI algorithm is used to obtain R component in each of two directions of upper left to lower right (labeled as [TeX:] $$\mathrm{D}_{1}$$ direction) and upper right to lower left (marked as [TeX:] $$\mathrm{D}_{2}$$ direction) at each [TeX:] $$\mathrm{L}_{\mathrm{B}}$$.

Step 2: The ARI algorithm is used at each of the [TeX:] $$\mathrm{L}_{\mathrm{G}}$$ directions in H and V directions to obtain R component.

Since the ARI algorithm is used along H and V directions in Step 2, which is the same as the interpolation method in Section 2.2, and then only the diagonal ARI algorithm of Step 1 is introduced.

Step 1 iteratively performs two kinds of algorithms (RI and MLRI algorithms) in two directions ([TeX:] $$\mathrm{D}_{1}$$ and [TeX:] $$\mathrm{D}_{2}$$ directions), and the RI algorithm flow of R component in [TeX:] $$\mathrm{D}_{1}$$ direction is shown in Fig. 2, where the dotted line direction represents [TeX:] $$\mathrm{D}_{1}$$ direction. The calculation process of MLRI algorithm in [TeX:] $$\mathrm{D}_{1}$$ direction and RI and MLRI algorithms in [TeX:] $$\mathrm{D}_{2}$$ direction are performed in the same way.

Firstly, an initial linear interpolation of R component in [TeX:] $$\mathrm{D}_{1}$$ direction is performed at [TeX:] $$\mathrm{L}_{\mathrm{B}}$$, and the interpolation results can be represented as:

where [TeX:] $$\tilde{R}_{(i, j)}$$ represents R component of image [TeX:] $$\tilde{R}$$ at [TeX:] $$(i, j)$$ , [TeX:] $$R_{(i-1, j-1)}$$ and [TeX:] $$R_{(i+1, j+1)}$$ represent R component of image R at [TeX:] $$(i-1, j-1)$$ and [TeX:] $$(i+1, j+1)$$ , respectively. Meanwhile, the superscript [TeX:] $$(\cdot)^{D_{1}}$$ represents that the interpolation process is in [TeX:] $$\mathrm{D}_{1}$$ direction, and the subscript [TeX:] $$\text { (.) }_{(i, j), 0}$$ represents the initial interpolation at [TeX:] $$(i, j)$$.

Secondly, at the kth iteration, an estimated image [TeX:] $$\breve{R}$$ is generated based on the linear transformation of the previous interpolation results, and R component of the image at [TeX:] $$(i, j)$$ can be expressed as:

where [TeX:] $$\left(a^{\mathrm{R}}, b^{\mathrm{R}}\right)$$ represents the linear coefficient of R component during the estimation process. Then, the residual component is obtained as:

where the superscript [TeX:] $$(\cdot)^{\mathrm{RD}}$$ represents the residual calculation process of R component in [TeX:] $$\mathrm{D}_{1}$$ direction. The residual value is linearly interpolated in [TeX:] $$\mathrm{D}_{1}$$ direction, and the result can be expressed as:

Finally, to obtain the results of the kth iteration interpolation in [TeX:] $$\mathrm{D}_{1}$$ direction, the residual of the linear interpolation is added to [TeX:] $$\breve{R}$$, which can be expressed as:

After the kth iteration, according to the steps of the ARI algorithm, the optimal iteration coefficients of the four algorithms need to be determined. Since there is only R component at this time, Eq. (2) can be rewritten as:

At this point, the optimal iteration coefficient [TeX:] $$k_{\text {best }}$$ obtained by Eq. (3) becomes:

where [TeX:] $$c_{(i, j), k}^{\mathrm{D}_{1}}$$ can be expressed as:

where the parameters [TeX:] $$(m, n)$$ are still set to [TeX:] $$(2,1)$$.

After obtaining the optimal iterative coefficients, [TeX:] $$\tilde{R}_{(i, j), k_{\text {best }}}^{\mathrm{D}_{1}}$$ needs to be adaptively combined. In this step, the direction interpolation results of the RI and MLRI algorithms are combined by weighted average:

where the superscripts [TeX:] $$(\cdot)^{D_{1}}$$ and [TeX:] $$(\cdot)^{\mathrm{D}_{2}}$$ represent the [TeX:] $$\mathrm{D}_{1}$$ direction and the [TeX:] $$\mathrm{D}_{2}$$ direction, respectively. The weights [TeX:] $$\omega_{(i, j)}^{\mathrm{D}_{1}, \mathrm{RI}}$$, [TeX:] $$\omega_{(i, j)}^{\mathrm{D}_{1}, \mathrm{MLRI}}$$, [TeX:] $$\omega_{(i, j)}^{\mathrm{D}_{2}, \mathrm{RI}}$$, and [TeX:] $$\omega_{(i, j)}^{\mathrm{D}_{2}, \mathrm{MLRI}}$$ in Eq. (12) can be expressed as:

4. Experimental Results and Performance Analysis

Obtaining high-quality restored images is the ultimate goal of demosaicking, so it is necessary to evaluate the quality of the interpolated images objectively as well as subjectively. This paper first introduces the test image and experimental environment, and then evaluates the image quality of the five algorithms, that is, RI, MLRI, IRI, ARI, and IARI, objectively and subjectively. Finally, the average interpolation calculation time of the five algorithms is utilized to reflect the complexity of the algorithm.

4.1 Test Image and Experimental Environment

To evaluate the quality of five algorithms, this paper selects two types test images from IMAX and Kodak datasets. The representative survey in [18] points out that these two datasets are reliable samples of the demosaick image quality test. The IMAX dataset contains 18 images of 500×500 pixels in TIF format, as shown in Fig. 3, while the Kodak dataset contains 24 PNG images with size of 768×512 pixels, as shown in Fig. 4.

The software used in all the simulation experiments in this paper is MATLAB R2016a, and the operating environment is a notebook computer with 2.40 GHz Intel Core i7-5500U CPU.

4.2 Comparison of Objective Quality of Images

In this section, the peak signal-to-noise ratio (PSNR) and the composite PSNR (CPSNR) are used to test the performance of the demosaicking algorithm.

As an important indicator of image quality, PSNR is denoted by [TeX:] $$Q_{\mathrm{PSNR}}$$ in this paper, which can be expressed as [11]:

where [TeX:] $$I_{\mathrm{in}}$$ is the input image and [TeX:] $$I_{\text {out }}$$ is the output image, and M X N is the size of the image. [TeX:] $$i$$ and [TeX:] $$ㅓ$$ represent the specific location of a pixel in the image. Since the PSNR is only for the evaluation quality of one of the R, G, and B colors, the overall evaluation of the image quality is not comprehensive. Therefore, this paper also introduces CPSNR, which is denoted by [TeX:] $$Q_{\text {CPSNR }}$$ in this paper, to evaluate the overall performance of the interpolation algorithm. The [TeX:] $$Q_{\text {CPSNR }}$$ can be expressed as [11]:

where represents the color plane of R, G, and B components, respectively.

Table 1 shows the comparison of the [TeX:] $$Q_{\mathrm{PSNR}}$$ of G, R, and B components after processing the 42 images in IMAX and Kodak datasets among five methods mentioned above.

From Table 1, it can be seen that the average [TeX:] $$Q_{\mathrm{PSNR}}$$ of each color component of the IARI algorithm is improved, compared with the other algorithms, in terms of the most images. Specifically, the average [TeX:] $$Q_{\mathrm{PSNR}}$$ of G component of the IARI and ARI algorithms is 1.25 dB, 0.81 dB, and 0.62 dB higher than that of the RI, MLRI, and IRI algorithms, respectively. In addition, the average [TeX:] $$Q_{\mathrm{PSNR}}$$ of R component of the IARI is 1.47 dB, 0.74 dB, 0.77 dB, and 0.12 dB higher than that of the RI, MLRI, IRI, and ARI algorithms, respectively. The average [TeX:] $$Q_{\mathrm{PSNR}}$$ of B component by the method of the IARI is 1.22 dB, 0.69 dB, 0.37 dB, and 0.29 dB higher than that of the RI, MLRI, IRI, and ARI algorithms, respectively. To improve readability, the graphs which is simulated by MATLAB R2016a are used to further illustrate, as shown in Figs. 5–7. These three graphs show the comparison of [TeX:] $$Q_{\mathrm{PSNR}}$$ of G, R, and B components among five algorithms, respectively. For convenience, Fig. 3(a)–(r) are numbered 1 to 18 in sequence, and Fig. 4(a)–(x) are numbered 19 to 42 in sequence. The image numbers are recorded as parameter α.

From Fig. 5, it can be seen that the curves of ARI and IARI algorithms overlap completely, that is, the IARI and ARI algorithms have the same effect on the interpolation of G component. From Figs. 6 and 7, it can be seen that in most cases, the curve of IARI algorithm is located above other curves, that is, for most images, the [TeX:] $$Q_{\mathrm{PSNR}}$$ of R and G components of the IARI algorithm are better than that of other four algorithms. Table 2 shows the comparison of [TeX:] $$Q_{\text {CPSNR }}$$, and the bolded data in table is the optimal results for each image. Moreover, to improve readability, Fig. 8 is used to illustrate the comparison of [TeX:] $$Q_{\text {CPSNR }}$$ among five algorithms.

Table 2 shows that the average [TeX:] $$Q_{\text {CPSNR }}$$ of the IARI algorithm is improved, compared with other algorithms, in terms of the most images. Specifically, the average [TeX:] $$Q_{\text {CPSNR }}$$ of the IARI algorithm is 1.43 dB, 0.68 dB, 0.42 dB, and 0.26 dB higher than that of RI, MLRI, IRI, and ARI algorithms, respectively. From Fig. 8, it can be seen that in most cases, the curve of IARI algorithm is located above other curves, that is, for most images, the [TeX:] $$Q_{\text {CPSNR }}$$ of the IARI algorithm is higher than that of other four algorithms.

In summary, experimental results show that the IARI algorithm is not all superior to the ARI algorithm in any case. This paper believes that the reason for these results is that in order to be consistent with the ARI algorithm of G component, the iteration parameters (m,n) in Eq. (11) are still set to (2,1) . However, (2,1) are not necessarily the optimal parameter in the interpolation process on R and B components. In summary, experimental results show that the interpolation effect of the IARI algorithm is better than that of other four comparison algorithms.

4.3 Comparison of Subjective Quality of Images

In this section, IMAX03 and IMAX08, as shown in Fig. 3(c) and (h), are selected to illustrate the subjective effects of the five algorithms, as shown in Figs. 9 and 10.

Comparing the flower diagrams of the original figure in Fig. 9(b) as a standard, in the complex marginal area where green grass and pink flowers cross each other, Fig. 9(c)–(e) grass and flowers have a certain degree of mosaic effect. The flower diagrams shown in Fig. 9(f) and (g) have a good recovery effect on the complex edge regions of this position, and the flower color recovery brightness of Fig. 9(g) is closest to Fig. 9(b). Comparing the wave point diagrams in the original figure of Fig. 10(b) as a standard, the recovery effects of the wave point areas of the red tie are not satisfactory in Fig. 10(c)–(e). Compared with the IARI wave point diagram shown in Fig. 10(g), the ARI wave point diagram of Fig. 10(f) basically restores the effect, but the black wave point of some positions produces a significant loss.

In summary, the subjective evaluation proves that the interpolation effect of the IARI algorithm is better than that of other four comparison algorithms.

4.4 Comparison of Algorithm Complexity

In the operating environment utilized in this paper, 42 images in IMAX and Kodak datasets are selected, and then each image is performed 15 times by using RI, MLRI, IRI, ARI, and IARI algorithms, respectively, and the average interpolation processing time is shown in Table 3. It can be seen from Table 3 that the average processing time of IARI and ARI is much longer than that of the other three algorithms, and the average time of IARI algorithm is 6.63 seconds longer than that of ARI algorithm. The experimental results show that the IARI algorithm has a higher computational complexity.

5. Conclusion

This paper focuses on ARI algorithm and illustrates the advantages and disadvantages of the four interpolation algorithms in image demosaicking by analyzing their principles, then proposes an improved algorithm IARI with better interpolation performance. According to the directional difference, IARI algorithm executes ARI algorithm twice on R and B components. By this way, the adaptive interpolation for R and B components is obtained. Experimental results show that the IARI algorithm achieves better performance than other four algorithms both in subjective assessment and objective assessment, especially in the complex edge area and color brightness recovery. The innovation of this paper is that the proposed IARI algorithm divides R and B components interpolation process into two categories according to the directional difference, then performs ARI algorithm interpolation twice according to different categories. Compared with the method by using ARI algorithm directly for both types of interpolation directions, the IARI algorithm is easier to operate. The improved algorithm proposed in this paper solves the problem of incomplete adaptive processing of R and B components in ARI algorithm interpolation process effectively. However, there are still some disadvantages in IARI algorithm. For example, the higher complexity of the IARI algorithm greatly sacrifices the interpolation speed of the algorithm, and the calculation time lasts a little long. Therefore, in the future, more research for improving the performance of IARI algorithm will be done. Under the premise of maintaining high performance, the complexity of IARI algorithm will be reduced as much as possible, therefore, the processing time of IARI algorithm will be reduced. In Eq. (11), the same parameters (m,n) = (2,1) as the calculated G component are selected, which might make the calculated values of R and B components not optimal and not accurate. Therefore, the setting of optimal parameters in R and B components processes should also be studied. The detailed future works are as follows. In [17], two important issues in the iterative RI process need to be further investigated: (1) the window size of the regression filter, and (2) the stopping criterion for the RI iterative process. Based on the reconstructed G component, the optimal parameters of R and B components should be determined using the description in [7], and then the optimal iteration coefficient will be determined.