Multispectral snapshot demosaicing via non-convex matrix completion

TL;DR

This paper introduces a non-convex matrix completion method for multispectral snapshot demosaicing, significantly improving image quality over existing techniques by accurately reconstructing missing spectral data.

Contribution

It proposes a novel non-convex matrix completion approach initialized with traditional algorithms to enhance multispectral image reconstruction from undersampled data.

Findings

Peak signal-to-noise ratio improved by 2-5 dB

Effective for p=16 mosaic sensors in diverse scenes

Outperforms current state-of-the-art methods

Abstract

Snapshot mosaic multispectral imagery acquires an undersampled data cube by acquiring a single spectral measurement per spatial pixel. Sensors which acquire $p$ frequencies, therefore, suffer from severe $1/p$ undersampling of the full data cube. We show that the missing entries can be accurately imputed using non-convex techniques from sparse approximation and matrix completion initialised with traditional demosaicing algorithms. In particular, we observe the peak signal-to-noise ratio can typically be improved by 2 to 5 dB over current state-of-the-art methods when simulating a $p=16$ mosaic sensor measuring both high and low altitude urban and rural scenes as well as ground-based scenes.

Tables

Table 1. Table 1 : Average PSNR over the 16 bands, with standard deviations. The best results for each image are highlighted in bold.

IMAGE	INIT		ASD	CGIHT	W2$\times$DCT
D0201	ID	$34.5\pm 0.8$	$36.3\pm 1.1$	$37.9\pm 1.1$	$34.3\pm 0.5$
	SD	$35.9\pm 0.7$	$38.1\pm 0.9$	$39.5\pm 1.2$	$35.8\pm 0.6$
D0301	ID	$36.0\pm 1.8$	$38.3\pm 2.2$	$39.2\pm 1.7$	$35.1\pm 1.1$
	SD	$37.3\pm 1.3$	$39.7\pm 1.4$	$40.5\pm 1.5$	$36.3\pm 1.0$
D0303	ID	$39.1\pm 1.0$	$41.2\pm 1.2$	$43.1\pm 1.4$	$38.8\pm 0.7$
	SD	$40.8\pm 0.8$	$43.4\pm 1.1$	$45.5\pm 1.3$	$40.5\pm 0.7$
D0307	ID	$34.1\pm 1.8$	$36.1\pm 2.3$	$36.9\pm 2.0$	$33.2\pm 1.1$
	SD	$35.2\pm 1.1$	$37.2\pm 1.2$	$37.8\pm 1.5$	$34.3\pm 0.8$
D2301	ID	$37.9\pm 0.3$	$39.4\pm 0.7$	$41.1\pm 1.3$	$37.8\pm 0.3$
	SD	$39.5\pm 0.3$	$41.5\pm 1.0$	$43.6\pm 1.5$	$39.7\pm 0.4$
AvLF	ID	$33.7\pm 0.7$	$35.6\pm 1.2$	$36.7\pm 2.3$	$33.5\pm 0.7$
	SD	$35.3\pm 0.7$	$37.8\pm 1.6$	$39.3\pm 2.6$	$34.7\pm 0.6$
StCh	ID	$36.2\pm 1.0$	$36.7\pm 1.7$	$36.6\pm 1.8$	$36.0\pm 1.0$
	SD	$37.3\pm 1.1$	$37.5\pm 1.0$	$37.2\pm 0.9$	$37.4\pm 1.0$
N04	ID	$37.3\pm 3.1$	$38.6\pm 2.6$	$38.2\pm 2.8$	$35.3\pm 1.9$
	SD	$35.6\pm 0.8$	$36.9\pm 0.8$	$35.4\pm 2.2$	$34.9\pm 0.9$
N06	ID	$33.3\pm 0.8$	$34.5\pm 0.6$	$35.2\pm 0.5$	$32.6\pm 0.4$
	SD	$33.5\pm 0.3$	$34.8\pm 0.4$	$35.3\pm 0.8$	$33.2\pm 0.2$
N08	ID	$33.5\pm 0.4$	$34.4\pm 0.6$	$35.6\pm 1.0$	$32.7\pm 0.2$
	SD	$33.2\pm 0.4$	$34.4\pm 0.7$	$35.5\pm 1.3$	$32.9\pm 0.3$
F05	ID	$36.1\pm 2.4$	$36.6\pm 1.9$	$36.5\pm 1.9$	$34.0\pm 1.4$
	SD	$35.0\pm 0.4$	$36.1\pm 0.8$	$35.7\pm 1.6$	$34.1\pm 0.7$
F06	ID	$40.0\pm 0.9$	$40.2\pm 0.8$	$39.9\pm 0.9$	$38.6\pm 0.7$
	SD	$38.9\pm 0.4$	$39.6\pm 0.9$	$39.3\pm 0.9$	$39.1\pm 0.5$
F07	ID	$34.9\pm 1.2$	$36.3\pm 1.5$	$36.1\pm 1.6$	$34.4\pm 0.9$
	SD	$34.6\pm 0.8$	$36.3\pm 1.2$	$35.1\pm 1.8$	$34.6\pm 0.8$
GB	ID	$34.3\pm 1.1$	$36.0\pm 0.9$	$36.7\pm 1.1$	$33.3\pm 0.5$
	SD	$34.9\pm 0.3$	$37.2\pm 0.9$	$37.5\pm 1.5$	$34.7\pm 0.7$
GD	ID	$34.1\pm 1.4$	$35.8\pm 1.1$	$36.2\pm 1.2$	$33.0\pm 0.5$
	SD	$34.6\pm 0.4$	$36.9\pm 0.9$	$36.6\pm 1.1$	$34.3\pm 0.7$
GR	ID	$35.3\pm 1.2$	$37.0\pm 1.0$	$37.3\pm 1.2$	$34.4\pm 0.5$
	SD	$36.0\pm 0.4$	$38.0\pm 0.9$	$38.1\pm 1.7$	$35.7\pm 0.6$

Equations

$$\underset{x}{\min }\Vert y-P_{\Omega }{\Psi }^{-1}x{\Vert }_2,\text{s.t.}\Vert x{\Vert }_0\le k,$$

$$\underset{X}{\min }\Vert y-P_{\Omega }X{\Vert }_2,\text{s.t.}\mathrm{rank}(X)\le r,$$

$${\text{PSNR}}_k=10{\log }_{10}\left(\frac{({\max }_{p\in \mathcal{P}}{I}_p^k{)}^2}{\frac{1}{|\mathcal{P}|}{\sum }_{p\in \mathcal{P}}(I_p^k-{\hat{I}}_p^k{)}^2}\right),$$

Full text

**Index Terms— ** Snapshot multispectral imaging. Sparse approximation. Compressed sensing. Matrix completion. Demosaicing.

1 Introduction

Multispectral imaging is the process of recording 2D arrays of information at multiple spectra (frequencies). Having access to such a rich, three-dimensional data cube allows different materials to be distinguished due to their differing spectral emission profiles. As a result, multispectral imaging is used in applications ranging from landmine detection, precision agriculture, and medical diagnosis to name but a few of its application domains; for a partial survey see the January 2014 special issue of IEEE Signal Processing Magazine [1]. The increased sensor size and acquisition time are some of the central obstacles to the more widespread use of multispectral imagery.

Snapshot mosaic multispectral sensors allow for a compact video rate multispectral imaging by acquiring only a fraction of the multispectral cube. For example111We illustrate the architecture through the IMEC sensor, but note there are numerous similar sensors such as the S 137 system by CUBERT., the IMEC SNm4x4 records 16 bands at a rate of 340 frames per second on a spatial two-dimensional $2048\times 1088$ pixel domain by only acquiring a single spectrum per pixel; specifically this is achieved by tiling the two-dimensional domain by $4\times 4$ pixel supercells where each supercell acquires the spectra independently.

Herein we demonstrate the efficacy of multiple methods for interpolating the missing values in the above snapshot mosaic data cube by simulating the undersampling from complete three-dimensional data cubes provided by DSTL as well as AVIRIS [2], Stanford SCIEN [3], Nascimento [4], Foster [5], IEEE GRSS Data Fusion Contest [6]. In addition to reviewing the existing state-of-the-art interpolation methods in Sec. 2.1, we demonstrate sparse approximation and matrix completion methods in Sec. 2.2 and 2.3 respectively, which we observe to substantially outperform the prior state-of-the-art. Over the above diverse data sets, we observe that non-convex compressed sensing and matrix completion methods initialised with traditional interpolations methods typically improve the peak signal-to-noise ratio by $2\text{\,}\mathrm{dB}5\text{\,}\mathrm{dB}$, see Table 1.

2 Algorithms for demosaicing

Demosaicing is the process by which the undersampled three-dimensional multispectral data cube has the missing entries approximated so as to simulate a full data acquisition. While most three-dimensional interpolations methods would be directly applicable, we consider a few methods previously used by the multispectral community, such as direct interpolation as described in Sec. 2.1 as well as sparse approximation regularisation methods in Sec. 2.2 and low-rank structure as presented in Sec. 2.3.

2.1 Direct interpolation methods

Brauers et al. [7] developed methods to estimate the missing values in the multispectral cube based on extending a spatial bilinear interpolation of the measured values to include any spectral correlation. The weighted bilinear interpolation (WB) for the $4\times 4$ pixel regular mosaic filter follows by padding the missing entries with zeros and convolving with the cartesian product of a discrete seven-pixel width filter $\frac{1}{4}[1\;2\;3\;4\;3\;2\;1]$. Then, the spectral correlation is included in the spectral difference (SD) method by a) taking the output of WB to independently compute, for each band, say $k$, the difference between the values of the measured pixels for spectrum $k$ and the WB interpolated values of every other band restricted to the support of the measured pixels of spectrum $k$, then b) applying WB to these spectral differences c) to form an approximation of the full spectrum $k$ by adding the output of step (b) to the difference with $l$ at the location of the measured pixels for spectrum $l$.

Mihoubi et al. [8] extended the SD approach to consider alternative ways to build correlations between the bands. In intensity difference (ID) they build spectral correlation by first constructing a spatial intensity map whose value at a pixel is the measurement for whichever spectra was measured at that spatial location, then this intensity map is averaged using a weighting based on the number of pixels per spectra contained in the averaging width. See Sec.3.2 for details on the choice of averaging used here and Mihoubi et al. [8] for a more general discussion. Hence, the difference between this averaged intensity map and the measurements for each spectrum is computed, the unknown values zero padded, and each band averaged such as in WB.

Interpolation methods designed in transform domains have been considered by Miao et al. [9] in the binary tree-based edge-sensing (BTES) method, which has the additional benefit of allowing for variable sampling densities per frequency band. However, we observe it is inferior to SD and ID described above in the setting of snapshot imaging. Pseudo-panchromatic image difference (PPID) [10] builds upon BTES and ID. However, due to the applicability of PPID to only some specific mosaic arrangements we leave comparison with our algorithms for a later time.

2.2 Sparse approximation and compressed sensing inpainting

Sparse approximation inpainting allows one to easily consider the interpolation of the under-complete snapshot data cube in transforms more general than the linear interpolation of (WB). In particular, one can assume that the image is well approximated by a sparse representation in a suitable image domain and exploit this structure to reconstruct it from undersampled measurements [11, 12, 13], e.g. by solving

[TABLE]

where $\Psi$ represents the transform in which the data is known to be compressible and $y$ is the full data cube projected by $P_{\Omega}$ to the undersampled locations. Degraux et al. [14] apply this model to a reconstruction of multispectral imagery acquired by mosaic snapshot cameras.

The primary challenge lies in two aspects: (i) the significant subsampling ratio of $1/K$, where $K$ is the number of spectral bands, and (ii) the selection of the suitable transform $\Psi$. The first problem can be overcome by initialising the state-of-the-art sparse approximation algorithms for solving (1) with sufficiently accurate initial estimates, such as those from the classical interpolation methods described in Sec. 2.1. As it was pointed out in [14], we find that, for natural scenes captured by snapshot multispectral imaging, a Kronecker product of 2D wavelet transform spatially and the discrete cosine transform for the spectral dimension is an effective choice for the representation $\Psi$. In particular, the 2D wavelet transform spatially includes elements of nearly global support to allow broad correlations as well as local elements to express fine detail and the discrete cosine transform models the slowly varying values in the spectral dimension.

2.3 Matrix completion

Rather than using local correlations, matrix (and tensor) completion exploit the correlation in the data cube through a low-rank structure, e.g. by solving

[TABLE]

where $y$ is the observed data, $X$ is a matrix corresponding to an unfolding of the complete three-dimensional data cube, and $P_{\Omega}$ is a restriction to the measured values as before.

Although the low-rank matrix completion problem is NP-hard in general (see [15] for a recent survey) there is a number of computationally fast solvers for the problem with provable convergence guarantees [16, 17, 18]. In fact, matrix completion has been previously applied to the reconstruction of subsampled multispectral imagery [19] by the Coded Aperture Snapshot Spectral Imager (CASSI) [20]. Here we show that matrix completion can be used also in the case of a more severe subsampling by snapshot mosaic camera designs if provided with suitable initialisation.

While there are many non-convex methods for matrix completion, we showcase two exemplary cases but expect that other non-convex methods would perform similarly well. We apply conjugate gradient iterative hard thresholding (CGIHT) [18] and alternating steepest descent (ASD) [17] to solve (2), providing them with an initial guess from either SD or ID. We show that both CGIHT and ASD can improve on the classical interpolation methods. This differs substantially from prior work both in terms of using more recent algorithms for matrix completion which have been shown to be more effective and initialising them with prior state-of-the-art interpolation methods, and moreover in that unlike [19] which treat each spectral band separately with $30\%$ undersampling, we vectorise the spatial dimensions to create a matrix of size $2,228,224$ by $16$ with $1/16$ undersampling. We observe that this unfolding which allows full correlation between the spectral information is particularly effective, often resulting in reconstructions which are visually indistinguishable from fully acquired data.

3 Numerical simulations

In this section, we show and explain the numerical results obtained by applying the methods discussed above.

3.1 Data sets

We consider the efficacy of the algorithms for demosaicing on the following data sets:

• •

High altitude airborne images from the AVIRIS [2] and 2018 IEEE GRSS Data Fusion contest [6]. AVIRIS line-scan captures $224$ spectral bands in the $380\text{\,}\mathrm{nm}2500\text{\,}\mathrm{nm}$ and the GRSS images have $48$ spectral bands in the range of $380\text{\,}\mathrm{nm}1050\text{\,}\mathrm{nm}$.

• •

Low altitute airborne images acquired at DSTL Porton Down, in August 2014, from which we selected $10$ representative radiance images of fields from a HySpex VNIR-1600 line-scan camera in the range $400\text{\,}\mathrm{nm}1000\text{\,}\mathrm{nm}$.

• •

Ground-based images from the Stanford SCIEN [3], Nascimento [4] and Foster [5]. The Stanford SCIEN images come from the line-scan HySpex VNIR-1600 camera.

We processed these data sets to simulate the spectrum measures by the IMEC SNm4x4 snapshot sensor with access to the complete data cube. Then, we undersampled the data cube following the sensor sampling pattern and the following simulations conducted.

3.2 Simulation setup

We implement and test recovery by two interpolation methods ID and SD, two matrix completion methods ASD and CGIHT and a compressed sensing version of CGIHT with a sparsifying transform as a Kronecker products of 2D Daubechies wavelets (W2) in the spatial and 1D Discrete Cosine Transform (DCT) in the spectral domain which we reference as W2$\times$DCT.

The iterative algorithms are terminated once the error in iteration $t$ is $\|P_{\Omega}X^{(t)}-y\|_{2}/\|y\|_{2}\leq 10^{-7}$ or at the $500^{th}$ iteration.

We report the quality of an image approximated by demosaicing by the peak signal-to-noise ratio (PSNR), defined as the log of the ratio between the maximum possible power of (the slide of) an image and the power of corrupting noise that affects the fidelity of its representation, computed in terms of the average squared difference (or mean squared error, MSE) between the reference image and its reconstruction:

[TABLE]

where $I^{k}$ and $\hat{I}^{k}$ are the $k$-th band slices of the reference cube and the reconstruction, respectively, and $\mathcal{P}$ denotes the set of all pixels.

We also employ the structural similarity (SSIM) index [21], which is a decimal value between $-1$ and $1$, with $1$ being reachable only in the case of two identical sets of data. SSIM is a perception-based model that considers image degradation as perceived change in structural information, while also incorporating important perceptual phenomena, including both luminance masking and contrast masking terms.

3.3 Results

Figure 1(a) shows the PSNR of each spectrum given its band centre, for a sample image from Nascimento [4], for the classical interpolation methods SD and ID as well as the compressed sensing (CS) and matrix completion techniques initialised with SD. Notice that, except for the first band, the matrix completion algorithms outperform SD and ID. On the other hand, the compressed sensing approach does not improve on the interpolation methods. In particular, note the overall incorrect contrast level resulting in yellowing of Figure 1(b). Moreover, we lose the sharpness of the edges in the balcony when employing CS W2$\times$DCT (Figure 1(b)), while we recover it with the matrix completion variant of CGIHT (Figure 1(c)).

To further emphasise how the different algorithms differ from each other we show the reconstructions PSNR from the corresponding reference images in Figure 2. Note in Figure 2(c) how ID accurately reconstructs the field image, taking into account the spectral correlation between the bands, but smooths the sharp details. The compressed sensing approach does a better job in identifying the edges (Figure 2(b)), but suffers from the same problem overall. On the other hand, the matrix completion CGIHT outperforms the other methods due to the presence of field-like uniformities.

Finally, as shown in Table 1, in the majority of cases the best performing algorithms are the matrix completion CGIHT and ASD initialised with SD, with improvements over both SD and ID from $2\text{\,}\mathrm{dB}5\text{\,}\mathrm{dB}$. StCh from Stanford SCIEN [3] seems to be the only outlier, with an improvement of just $0.2\text{\,}\mathrm{dB}$.

Being directly related to the number of bands, the rank of the spectral unfolding seems to be effective in capturing the spectral information of the analysed datasets. Our results suggest a high correlation between the frequency bands and a low-rank structure of the spectral unfolding of our images, in which most of the information is contained in the first 3 singular values of the spectral unfolding.

4 Conclusions

We provide a numerical comparison of multispectral demosaicing by traditional interpolation, sparse approximation and matrix completion methods. Our experiments demonstrate that non-convex matrix completion typically improves reconstruction by $2\text{\,}\mathrm{dB}5\text{\,}\mathrm{dB}$ over the current state-of-the-art methods. This differs substantially from prior work in terms of employing matrix completion on the spectral unfolding of the image in the context of demosaicing, initialising it with classical interpolation methods and using more recent non-convex matrix completion algorithms.