Generic Multispectral Demosaicking Based on Directional Interpolation

The low-cost snapshot multispectral spectral imaging systems with multispectral filter array (MSFA) require generic MSFA demosaicking methods to generate the multispectral image (MSI) with the variable number of spectral bands depending on the applications. Most of the existing MSFA demosaicking methods are either non-generic or perform inadequately. This paper presents a new generic MSFA demosaicking method based on the directional weighted interpolation. Our proposed method calculates the four directional estimates around the location of the unknown pixel and combines them in a weighted manner using the local edge magnitude in the corresponding directions to estimate the missing pixel values. Experimental results confirm that the proposed demosaicking method provides improvement, compared to the state-of-the-art generic demosaicking methods in terms of both subjective and objective evaluations on the two benchmark MSI datasets.


I. INTRODUCTION
M ULTISPECTRAL images (MSIs) are rich in the spectral and spatial informational about the scene than standard RGB images. Hence, MSIs are more useful than RGB images in the area of remote sensing [1]- [3], satellite imaging [4], medical imaging [5]- [7], food industry [8]- [10], and computer vision [11]- [17]. MSIs are generally acquired by multispectral spectral imaging systems, which use either multiple exposures with single imaging sensor or single exposure with multiple imaging sensors. Use of multiple exposures make multispectral imaging systems unsuitable for moving objects and use of multiple senors make these systems of high cost. Motivated by the use of a color filter array (CFA) to capture RGB images, a multispectral filter array (MSFA) based imaging system simultaneously obtains all spectral bands information from a single shot using a single imaging sensor. Its capability to capture MSI using a single senor opens up a vast range of possibilities in various applications.
The image taken by a single sensor based multispectral imaging system stores just one band's information at each pixel location. This captured raw image using an MSFA is called an MSFA image, similar to the CFA image in the RGB domain. An interpolation method called MSFA demo-saicking is used to generate the full MSI from the MSFA image similar to the CFA demosaicking in the color (RGB) domain. The quality of the generated MSI not only relies on the MSFA demosaicking algorithm but also on the MSFA patterns used to capture the initial raw image. There are many generic MSFA patterns (non-redundant MSFA pattern [18], [19], uniform MSFA pattern [20], binary tree based MSFA pattern [21], and random [20]) are given in the literature, but binary tree based MSFA patterns are most preferred in the multispectral domain due to their compact nature as compared to other MSFA patterns [22], [23].
Several efficient CFA demosaicking methods [24], [25] have been proposed, but their direct extension to MSFA demosaicking is neither possible nor straightforward. MSFA demosaicking is difficult because the spectral and spatial correlation properties of MSIs differ greatly from those of RGB images. The spectral bands in the MSFA image are sparsely sampled, particularly in the MSFA image with a higher number of bands. Because of the sparsity, the MSFA image has few pixels of the same band, which lessens the spatial correlation. The weak spectral correlation of the band with the distant bands also makes MSFA demosaicking even more challenging [26]. The current MSFA demosaicking methods can be categorized into two classes. The first class contains the MSFA demosaicking methods that are limited to the specific number of bands MSIs [23], [27]- [30] and the second class contains generic MSFA demosaicking methods [18]- [20], [22], [26], [31]- [33] which can be used to generate any K-band MSI from the MSFA image. The generic demosaicking methods are picked due to their broad applicability and ability to provide multispectral imaging system manufacturers with numerous customizing choices through a single procedure.
In this paper, we propose a generic MSFA demosaicking method that generates the complete MSI from the MSFA image captured using binary tree based MSFA patterns. Several other existing generic MSFA demosaicking methods [22], [31]- [33] use MSFA images generated using binary tree based MSFA patterns. All these methods, except [31], use both the spatial and spectral correlations and their efficacy depends on the better utilization of spectral correlation present in the MSFA image. We aim to better utilize the spectral correlation and considering the insights presented in [34] (constrained to Bayer CFA only), a new MSFA demosaicking method is proposed that uses the multiple directional estimates of bands based on the spectral correlation to calculate the missing pixel values. In the case of color (RGB) images, bands have higher PoA and fixed pixel arrangements around the unknown pixel locations. However, in the multispectral domain, bands have lower (each band has PoA equal to 1 16 in the considered 16-band MSFA image) and varying PoA depending on the number of bands in the MSFA image. Different pixel arrangements of bands w.r.t. each other, are also possible depending on the number of bands in the MSFA image and the PoA of each spectral band. So directly applying [34] demosaicking method to the multispectral domain is not possible due to the much lower PoA of bands in the MSIs, and it is not straightforward to propose a generic MSFA demosaicking method in the multispectral domain due to the reasons mentioned above. In work [35], authors extended [34] to the multispectral domain; however, only for the specific band size MSIs, where number of bands is 2 p × 2 p (p = 1, 2, ...), as the non-redundant MSFA patterns and binary tree based MSFA patterns are same in these cases. Otherwise, their approach is not applicable to the compact and preferred binary tree based MSFAs. Further, [35] uses an existing method for the initial estimation of the multispectral image, so the steps of [34] were easily extendable. But the use of an initialization method also leads to bias and has the cascading effects of error. On the contrary, the proposed method uses the intelligent interpolation scheme based on the PoA of spectral bands without relying on any other method for an initial estimation.
The proposed approach uses the binary tree based MSFAs and spectral correlations effectively and the experimental results on two standard multispectral datasets illustrate its better performance compared to existing generic demosaicking approaches. The major contributions of the proposed MSFA demosaicking method are as follows: (i) The proposed demosaicking method is generic; therefore, it can be used to generate MSI from any K-band MSFA image captured using binary tree based MSFA pattern. (ii) To interpolate missing pixel values of a band, the proposed demosaicking method calculates the multiple directional estimates of that band based on spectral correlation and combines these estimates in a weighted way based on the local edge magnitude in the corresponding direction. (iii) The proposed demosaicking method uses a specific order of pixel locations of bands to interpolate each band at the locations of other bands based on the PoA of bands and the binary tree used to generate the MSFA patterns. (iv) The proposed method is progressive and chooses the interpolation scheme intelligently based on the PoA of bands, and the proposed method does not depend on any other existing method for the initial estimation of MSI.
The rest of the paper is organized as follows: Section II briefly discusses the existing MSFA demosaicking methods. The proposed demosaicking method is described in section III. Experimental results are shown in Section IV, and Section V concludes the paper.

A. BAND-SPECIFIC MSFA DEMOSAICKING METHODS
Monno et al. [27] proposed a 5-band MSFA demosaicking method. This method used binary tree based MSFA patterns to capture the 5-band MSFA image. In this 5-band MSFA patterns, the authors kept the PoA of G band equals to 1 2 . To generate the full MSI from the MSFA image, the authors first estimated the G band employing the adaptive kernel and later used the interpolated G as the guide image to estimate other bands. [29] used the MSFA pattern given by [27] and proposed an MSFA demosaicking method for 5-band MSIs based on frequency domain analysis. Mihoubi et al. [36] utilized the idea of intensity image, defined as the average of all the bands and considered strongly correlated with each band than bands taken pairwise. Their proposed method is non-generic and constrains to square-shaped non-redundant MSFA patterns. Later, in work [28], Mihoubi et al. presented a 16-band MSFA demosaicking method based on the concept of pseudo panchromatic image. Sun et al. [23] proposed an MSFA demosaicking method to generate the full MSI for the 9-band MSFA image captured using binary tree based MSFA pattern. To generate the MSI, the authors first interpolated the middle band (fifth band) using the neighboring known pixels value to the unknown pixel location using the image gradients and later used the interpolated fifth band as the guide filter to interpolate other bands.

B. GENERIC MSFA DEMOSAICKING METHODS
Miao et al. presented a generic binary tree based MSFA patterns [21] and the binary tree based edge sensing (BTES) generic MSFA demosaicking method [31]. BTES method employs only spatial correlation, and the spatial correlation decreases as the number of spectral bands increases in the MSFA image. Therefore, BTES shows poor performance on the higher band images. Brauers and Aach [18] proposed weighted bilinear (WB) and spectral difference (SD) interpolation based MSFA demosaicking method to generate 6band MSI from the MSFA image capture using 6-band nonredundant MSFA pattern arranged in 2 × 3 matrix. Later, WB and SD are generalized to generate any K-band MSIs using non-redundant MSFA patterns by [19] and [37], respectively. Further, Mizutani et al. [26] proposed an iterative spectral difference (ISD) based MSFA demosaicking method (improved version of [18]) by iteratively repeating the interpolation process. The number of repetitions is determined by the correlation between spectral bands. The number of repetitions in strongly correlated spectral bands is greater than in loosely correlated spectral bands.
In [32], the authors presented a generic MSFA demosaicking method (PBSD) for the binary tree based MSFA images by progressively using WB and SD. In work [22], the authors formulated new convolution filters based on the PoA of bands in binary tree based MSFA patterns and presented a generic MSFA demosaicking method (PCBSD) that utilized these PoA based filters to perform weighted interpolation and spectral difference based interpolation to calculate the missing pixel values. Recently, Gupta et al. proposed a generic MSFA demosaicking algorithm (APMID) [33] utilizing the binary tree based MSFA pattern. In the considered MSFA patterns, the authors kept the PoA of the middle band higher than other bands as the middle band is used to evaluate all other bands. In APMID, the authors first estimated the middle using adaptive and progressive interpolation inspired by the adaptive color image demosaicking [38]. Later, other bands are interpolated using the progressive spectral difference method by utilizing the interpolated middle band. On the higher band MSIs, APMID does not perform satisfactorily as its performance relies upon the middle band interpolation. The PoA of the middle band becomes lower in a higher number of bands MSIs, and the APMID utilizes only spatial correlation to interpolate the middle band at some locations in the higher bands MSIs.
Few deep learning-based MSFA demosaicking methods [39]- [41] are recently proposed. Some of these methods [39], [40] are constrained to specific band-size images and may require a complete architectural change and, therefore, require to re-train to make these methods applicable to different band-size MSIs. Even these deep learning based methods need some of the original MSIs to train and validate the model, and the original images are not accessible in the real-time imaging scenario.

C. OUR CONSIDERED BINARY TREE BASED MSFA PATTERNS
MSFA patterns decides the basis of the raw information captured in the MSFA image, and they are critical to the effectiveness of MSFA demosaicking methods. In this paper, we utilize binary tree based MSFA patterns given by [21]. Binary tree based MSFA patterns can be designed to capture any band-size MSFA image. Many MSFA demosaicking methods [27], [29], [42]- [47] effectively use these MSFA patterns. We show the binary tree based MSFA pattern used in the proposed method to capture the MSFA image of 5-16 band in Fig. 1.

III. PROPOSED MSFA DEMOSAICKING METHOD
This paper proposes an MSFA demosaicking method that takes advantage of both spectral and spatial correlation present in the MSFA image. MSFA image is captured using a single sensor with binary tree based MSFA and has sparse information of the spectral bands. The proposed demosaicking method is generic; therefore, it can be applied to generate any K-band MSI from the respective MSFA image. Our proposed method exploits the properties of the binary tree based MSFA patterns where all the bands with similar PoA, irrespective of the different number of bands in the MSIs, have similar pixel arrangements. Our proposed method uses the PoA of bands in the binary tree based MSFA patterns to choose the interpolation scheme intelligently. The proposed method calculates four directional estimates of the band around the unknown location and combines them based on the weights in respective directions based on the interpolation scheme. These weights are computed by utilizing the edge magnitude information from the band chosen for interpolation in the respective direction, whereas [34] uses both bands (the band (k) that is chosen for interpolation and the band (l) at which location band k is to be interpolated) in the weight calculation and band l may not be available due to its lower PoA compared to the band k in the multispectral domain.
The proposed method is progressive, i.e., it first estimates the part of missing pixel values. Later, it uses the estimated values with the initially known pixel values in the MSFA image to estimate other unknown pixel values. The proposed method is independent of the order of the bands taken for interpolation. However, for any band selected for the interpolation process, the ordering of pixel locations picked for the interpolation is fixed and different for each band. This order relies upon the binary tree utilized to generate the MSFA pattern and on the PoA of the bands. Due to the progressive nature of the proposed method, the PoA of each band is updated after interpolating the missing pixel values at other band pixel locations, and the proposed method chooses the correct interpolation scheme based on its updated PoA to interpolate at pixel locations of remaining bands, as discussed later in the following subsections. Algorithm 1 illustrates the steps of the proposed MSFA demosaicking   method. The proposed MSFA demosaicking method mainly has the following two components.

A. PIXEL SELECTION SCHEME
To interpolate any band k, our proposed method follows a particular ordering of unknown pixel locations selected for the interpolation process. This order of selection of unknown pixel locations picked for interpolation of any band k is given by the pixel selection scheme corresponding to band k.
The pixel selection scheme utilizes the binary tree employed to form the MSFA pattern and starts from the leaf corresponding to the band selected for the interpolation in that binary tree. The pixel selection scheme first interpolates at the pixel locations of the sibling band of the band selected for the interpolation, in the binary tree. Then the scheme moves one level in the tree and selects the pixel locations corresponding to the sibling band of its parent. This pixel selection scheme progresses till the root of the binary tree is visited. If any sibling band is an internal node at any stage during interpolation, the bands corresponding to the leaves under that sibling node are considered irrespective of their order.
For example, we take 5-band MSFA image as shown in Fig. 2(a) and binary tree used to construct 5-band MSFA pattern shown in Fig. 2 shown in Fig. 2(e) and now PoA of band 1 becomes 1 as shown in Fig. 2(i).

B. INTERPOLATION SCHEME
The missing pixel values of band k at the location of band l at (i, j) is estimated using the spectral correlation of band k and band l in the neighborhood of (i, j). The interpolation scheme used to interpolate any band k depends on the PoA of band k. Here, we consider MSIs from 5-16 band, the initial PoA of any band k (P oA(k)) will be from the set { 1 2 , 1 4 , 1 8 , 1 16 } as per our considered MSFA patterns. As the proposed method is progressive, the PoA of band k is updated during the interpolation process as defined in the pixel selection scheme. Based on the PoA of band k at any stage of the interpolation scheme, we define four possible cases of the pixel arrangements of band k w.r.t unknown center pixel location (at the location of band l) in our proposed method as shown in Fig. 3. The interpolation scheme of band k depends on these four cases as described in the following section of the paper.

1) Case A and Case C
These cases A and C are applicable when P oA(k) is greater than equal to 1 2 and between [ 1 8 , 1 4 ), respectively during any stage of the interpolation scheme. The missing pixel values of band k at the center location (i, j) is estimated along four directions (north (N), south (S), east (E), and west (W)) as shown in Figs 3(a) and 3(c) by utilizing the spectral correlation between band k and band l. We consider t = 1 for Case A and t = 2 for Case C in the following equations.
where I k i,j is pixel value of band k at pixel location (i, j). We use weighted interpolation method to combine these four directional estimates. The weights along these four directions are calculated using edge magnitude of band k in these directions. So the larger the edge magnitude, the smaller the contribution from that directional estimate. The weight along these directions are estimated as follow: Finally, band k value is estimated as weighted sum as follow: Similarly, we calculate the weights along these four diagonal directions based on the edge magnitude of band k in these VOLUME x, 2022 Compute corresponding weights in these directions using Eqs. 5-8.

13
Estimate the value ofĨ k using Eq. 18. 14 end directions. The weight along these directions are estimated as follow: Finally, band k value is estimated as weighted sum as follow: Algorithm 1 summarizes the proposed MSFA demosaicking algorithm. Algorithm 2 summarizes the interpolation scheme for band k at the pixel locations of the band l and In Algorithm 1, it may be noted that to interpolate band k at the pixel locations of band l (selected from the pixel selection scheme for band k), P oA(k) should be less than four times the PoA of band l (P oA(l)). Otherwise, due to the small PoA of band l compared to the PoA of band k, the required number of known values of band l will not be available in the local neighborhood considered to estimate the band k at the pixel locations of band l, as per the cases mentioned in the interpolation scheme. In such cases, Algorithm 1 breaks the interpolation procedure for band k for the time being, and it proceeds ahead for interpolating other bands (including band l), leading to the enhancement of PoAs of these other bands. As the k and likely few other bands are not completely estimated in the first major iteration, Algorithm 1 continues the next iteration and would be reconsidering interpolating band k at the pixel locations of band l. As the P oA(l) has been increased, band l will be sufficiently available now to support the estimation of band k at the locations of band l. This process continues iteratively for all bands and finally, Algorithm 1 generates the estimated MSIĨ.

IV. EXPERIMENTAL RESULTS
We compare the performance of our proposed MSFA demosaicking method with existing generic MSFA demosaicking methods on MSIs from the CAVE [48] and TokyoTech [27] datasets. The CAVE dataset has 31-band MSIs of size 512 × 512. These MSIs are captured from 400nm to 700nm at every 10nm spectral gap. The TokyoTech dataset also contains 31-band MSIs of size 500 × 500 and these MSIs are capture from 420nm to 720nm at every 10nm spectral gap. We compare our proposed MSFA demosaicking method with WB [19], SD [18], BTES [31], ISD [26], PBSD [32], PCBSD [22], and APMID [33]. We use PSNR, sRGB PSNR, and CIEDE2000 [49] as image quality metrics for the objec- This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication.   [18], [19], SD [18], BTES [31], ISD [26], PBSD [32], PCBSD [22], and APMID [33]) on images from the TokyoTech dataset.  tive evaluation. For the subjective evaluation, we transform the MSI into the sRGB image [27] and compare different MSFA demosaicking methods based on the visual artifacts present in the sRGB images. We compare the performance of different generic MSFA demosaicking methods on the MSIs with the number of bands (K) ranging from 5 to 16 as many existing manufacturers [50]- [52] develop multispectral imaging systems that capture MSIs generally from 5-band to 16-band. In our experiments, to have the K-band reference image from the 31-band image, we select K bands at uniform spectral gaps starting from the 1 st band of the 31-band image. First, we use mosaicking to get the MSFA image, and later, we perform demosaicking to reconstruct the MSI. Finally, the reconstructed MSI is compared to the corresponding reference MSI in order to assess the effectiveness of the MSFA demosaicked methods.
Tables 1 and 2 show the PNSR values of different MSFA demosaicking approaches on the MSIs of the CAVE and TokyoTech datasets, respectively. Our proposed method consistently surpasses other generic demosaicking approaches almost by 1 dB on the lower number of bands MSIs and by 2 dB on the higher number of bands MSIs on both the datasets. Our proposed method exceeds the second best MSFA demo-  [18], [19], SD [18], BTES [31], ISD [26], PBSD [32], PCBSD [22], and APMID [33]) on images from the TokyoTech dataset.   [18], [19], SD [18], BTES [31], ISD [26], PBSD [32], PCBSD [22], and APMID [33]) on images from the CAVE dataset. saicking method (APMID) by 1.46 dB and 2.20 dB on the CAVE and TokyoTech datasets, respectively as average over 5-band to 16-band MSIs. All MSFA demosaicking methods (BTES, PBSD, PCBSD, APMID, and Ours) based on binary tree based MSFA patterns outperform non-redundant MSFA patterns based demosaicking methods (WB, SD, and ISD). This is due to the compact nature of binary tree based MSFA patterns as compared to the non-redundant MSFA patterns [22]. Even, the performance of the WB, SD, and ISD methods significantly degrade on the prime number (e.g., 5,7,11,13) of bands MSIs as compared to the non-prime number (e.g., 6,8,9,12) of bands MSIs due to the highly noncompact nature of prime number of bands MSFA pattern as compared to the non-prime number of bands MSFA patterns. Table 7 shows the performance of different MSFA demosaicking methods on the individual MSIs of 10-band and 15-band from the CAVE dataset. Clearly, our proposed method generates almost all MSIs with better quality from the corresponding MSFA images than other generic MSFA demosaicking methods. Similarly, Table 8 [18], [19], SD [18], BTES [31], ISD [26], PBSD [32], PCBSD [22], and APMID [33]) on images from the TokyoTech dataset. all MSIs from the TokyoTech dataset. Fig. 4 shows the performance of individual band of different band-size MSIs form both datasets. Our proposed method outperforms other demosaicking methods on all bands on 8-band MSIs and on the all bands (expect 9 th ) on 14-band MSIs from the TokyoTech dataset as shown in Fig. 4 (c,d). Our proposed method also performs better on all bands (except first band) on 10-band and 15-band MSIs from the CAVE dataset than other methods as shown in Fig. 4 (a,b). We transform the K-band images to the sRGB domain to assess colorimetric correctness. Tables 3 and 4 show the PSNR values of sRGB images converted from MSIs generated by the different MSFA demosaicking methods. Our proposed MSFA demosaicking method consistently outperforms other generic demosaicking methods by more than 1 dB on the sRGB images converted from all band-size MSIs on both the datasets. On average over all band-size images, our proposed method exceeds the second best MSFA demosaicking method (PCBSD) by 1.35 dB and 2.07 dB on the CAVE and TokyoTech datasets, respectively. Our proposed MSFA demosaicking method performs better than other demosaicking methods in terms of CIEDE2000 metric, as shown in Tables 5 and 6 on the CAVE and TokyoTech   8   VOLUME x, 2022 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and   [18], [19], SD [18], BTES [31], ISD [26], PBSD [32], PCBSD [22], and APMID [33]) on individual images of the CAVE datasets. Images  WB  SD  BTES  ISD  PBSD  PCBSD APMID  Ours  WB  SD  BTES  ISD  PBSD  PCBSD   dataset, respectively. We also examine the performance of the proposed method by considering the weights used in intermediate steps of [35] and the comparative results are presented in Table 9. As weights defined in [35] use both the bands (the band k whose value to be interpolated and band l at which location band k to be interpolated) and band l values may not be available due to its lower PoA than PoA of band K, we also use the estimated values of band l using the WB method similar to [35]. Clearly, our proposed method with our defined weight performs better than by using weights of [35] [18], [19], SD [18], BTES [31], ISD [26], PBSD [32], PCBSD [22], and APMID [33]) on images from the TokyoTech dataset.    [18], [19], SD [18], BTES [31], ISD [26], PBSD [32], PCBSD [22], and APMID [33]) on different band size MSIs of both datasets.  Table 10 represents the average running of different MSFA demosaicking methods on MSIs from the CAVE and Toky-oTech datasets. We use MATLAB 2018b to implement these methods and execute them on a 64-bit MSI machine with a Core-i7-8750H processor with 8 GB random access memory and 2.20 GHz processing speed with Windows 10 Home as the operating system. WB and APMID are the fastest methods compared to other methods. However, our proposed method achieves the best results on image quality metrics than all other demosaicking methods considered and it gives better or nearly the same computation times than that of other demosaicking methods like ISD, PBSD, SD, and PCBSD on both the datasets. performance on the Cloth, Party, and Character images, respectively from the TokyoTech dataset and Figs 7 and 8 indicate performance on the Cloth and Egyptian Statue images from the CAVE dataset. Both WB and BTES methods generate blurred images, and blurriness increases as the number of bands increases in MSIs. PBSD, PCBSD, and APMID produce significant artifacts around the edges. Our proposed method generates good quality images with lesser artifacts than other demosaicking methods. VOLUME x, 2022 11 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and

V. CONCLUSION
In this paper, we present a generic MSFA demosaicking method that can reconstruct the MSI from the highly sparse MSFA image. The MSFA images are captured using binary tree based MSFA patterns highly preferred in the MSFA demosaicking. The proposed method does not reply on any other method for intial estimate of the MSI. To estimate the missing pixel values of any band, we present a progressive interpolation scheme that uses the PoA of that band during the interpolation process to select an appropriate case for interpolation. Using that appropriate case, we calculate four directional estimates of that band around the unknown location and combine these four estimates based on weights calculated using the edge magnitude in those respective directions. The experimental results confirm that our proposed MSFA demosaicking method outperforms other generic MSFA demosaicking methods both quantitatively and qualitatively on two benchmark datasets.