Compressive Hyperspectral Imaging with Side Information

I Introduction

Hyperspectral imaging techniques have been widely applied in various fields such as astronomy , remote sensing , and biomedical imaging . Unlike ordinary imaging systems, hyperspectral imagers capture a three-dimensional (3D) datacube, i.e., a 2D array of vectors that contain the spectral information at each spatial location. Inspired by compressive sensing (CS) , researchers have adopted joint sensing and reconstruction paradigms that measure a subset of the 3D spectral datacube and utilize CS inversion algorithms to retrieve a 3D estimate of the underlying hyperspectral images. This sensing paradigm follows the traditional benefits of CS - reduced data rates and system complexity at the expense of computational algorithmic development and postprocessing.

The coded aperture snapshot spectral imaging (CASSI) systems are examples of CS hyperspectral imaging systems. The CASSI paradigm encodes each of the datacube’s spectral channels with a unique 2D pattern, which is the underlying operating principle behind code division multiple access (CDMA). CASSI systems form an image onto a coded aperture placed at an intermediate image plane to spatially modulate the datacube with high-frequency patterns (Fig. 1). A disperser placed in a pupil plane behind the coded aperture spectrally shifts the coded image, effectively granting each spectral channel onto its own unique coding pattern when multiplexed onto a monochrome sensor. Unlike other hyperspectral imaging techniques , CASSI can obtain a discrete 3D estimate of the target spectral datacube from as little as a single 2D measurement. Such operation renders the system’s forward model highly underdetermined; inversion requires use of CS algorithms . Compressive hyperspectral imaging requires the signal under estimation to be sparse in a basis that is incoherent with the system sensing matrix . The reconstruction is accomplished using optimization algorithms, such as gradient projection for sparse reconstruction (GPSR) or two-step iterative shrinkage/thresholding (TwIST) . GPSR assumes sparsity of the entire datacube in a fixed (wavelet) basis, while TwIST is based on a piecewise-constant spatial intensity model (when the total variation norm is used) for hyperspectral images.

Distinct from the above optimization algorithms, blind CS algorithms have been applied to CASSI systems by learning dictionaries from the measurements. Blind CS inversion strategies seek to recover 3D patches of the hyperspectral datacube jointly with the shared dictionary (inferred from the measurements). Each patch is a sparse combination of the dictionary atoms. Since the dictionary is unknown a priori, this is called blind CS . This blind CS model shares statistical strengths among similar image patches at different spatial locations. Additionally, this CS approach is task-driven since the learned dictionary is appropriate for different tasks . In this paper, we propose a new blind CS model that imposes compressibility , rather than sparsity, on the recovered dictionary coefficients. Specifically, we use global-local shrinkage priors on the dictionary coefficients for each patch, under a Bayesian dictionary learning framework. This setting drives the reconstructed dictionary coefficient vector to contain many small (i.e., close to zero) components while allowing for a few large components, but avoiding explicit sparsity enforcement. This model is feasible for high-dimensional hyperspectral signals since it can extract more information from the limited (in situ learned) dictionary atoms.

The quality of the reconstructed hyperspectral images relies on the conditioning of the sensing matrix. More accurate recovery is possible with additional measurements of the same scene or by using digital micromirror (DMD) arrays ; however, these methods increase system cost and energy consumption. This paper features a blind CS algorithm that employs the RGB image of the target scene (as side information) and demonstrates substantial improvements on reconstruction quality with a single CASSI measurement. Particularly, the proposed algorithm learns a joint dictionary from a single CASSI measurement and the corresponding RGB image (of the same scene) and then reconstructs the hyperspectral datacube. Since the RGB image is strongly correlated with the hyperspectral images, this joint model will dramatically aid the dictionary learning algorithm, thereby improving the reconstruction.

Furthermore, we propose a new camera using spatial-light modulation (SLM) for an active coding compressive spectral imager, to jointly modulate the spatial and spectral information of the scene on the intermediate image plane. This technology differs from CASSI and DMD-modulated CS imagers in that no dispersive element is required for multiplexing the spectral channels.

The reminder of the paper is organized as follows. Section II reviews the mathematical model of CASSI and introduces the new camera. Section III develops the new blind CS algorithm to reconstruct hyperspectral images. Experimental results are presented in Section IV, and Section V summarizes the paper.

II Hardware

In this section, we first review the CASSI imaging process and then introduce the new camera. Finally, a shared mathematical model for both cameras is proposed.

A common CASSI design in Fig. 1 adapts the available hardware to encode and multiplex the 3D spatiospectral information $f(x,y;\lambda)$ onto a 2D detector. This modulation process is based on physically shifting the hyperspectral datacube multiplied by the binary coding pattern $T(x,y)$ via chromatic dispersion. The coding pattern consists of an array of square, binary (100% or 0% transmission) apertures, which provide full or zero photon transmission. The camera’s disperser (in traditional hyperspectral CS , this is a prism or diffraction grating) laterally displaces the image as a function of wavelength $d(\lambda-{\lambda}_{c})$ of the coded aperture pattern, where $d(\lambda)$ is the lateral wavelength-dependent shift and ${\lambda}_{c}$ represents the disperser’s central wavelength. As previously mentioned, this coding process can be considered a form of CDMA, whereby each channel is modulated by an independent coding pattern on the detector plane. The detector integrates the spectrally-shifted planes along the spectral axis. Information can be recovered by separating channels based on their projected coding patterns. This multiplexing process can be represented as the following mathematical forward model:

where $g(x,y)$ is the continuous form of the detector measurement, representing the sum of dispersed coded images. Since the detector array has been spatially pixelated by the detector pixel pitch $\Delta$, the discreatized detector measurement becomes:

Finally, the discrete measurement for each pixel can be illustrated as:

where $T_{m,n}$ is the spatially encoding pattern, $f_{m,n,j}$ represents the discretized spectral density of $f(x,y;\lambda)$, and $\epsilon_{m,n}$ is the noise.

II-B A New Camera: SLM-CASSI

Here we use an SLM to encode the 3D spatial-spectral information on a 2D gray-scale detector - a similar CDMA strategy to CASSI for snapshot hyperspectral imaging. An SLM is an array of micro cells placed on a reflective layer; each cell has a nematic liquid crystal that adds pixel-level phase to the incident light as a function of voltage, thereby changing the polarization state on a per-pixel basis . Each layer of the LC can be considered as a thin optical birefringence material; the orientation and the relative refractive index difference between the fast and slow axes determines the effective birefringence. Since most birefringent phase modulators are nominally sensitive to wavelength, this element can assign wavelength-dependent transmission patterns to multiplex every spectral channels for the compressive measurement. The hyperspectral slices can be separated from the coded data via CS inversion algorithms.

A schematic of this SLM-based CASSI system is shown in Fig. 2. An objective lens (L1) images the scene from a remote distance, and then the collimation lens (L2) and imaging lens (L3) relay the image through the polarizing beamsplitter and the achromatic quarter-wave retarder onto the SLM. The retarder serves to increase the contrast of the SLM polarization array by compensating the extra phase retardation in the liquid crystal device. The spatiospectral coding function is implemented on the SLM as an 8-bit, pseudo random, pixel-level grayscale phase retardation pattern. The phase retardation manifests as wavelength-dependent ampitude modulation upon re-entry through the polarizing beamsplitter toward the camera. Finally, the modulated image is projected by the imaging lens (L4) onto the detector plane and then recorded by the detector. The following mathematical forward model can be used to describe the multiplexing process of the compressive sampling:

where $T(x,y,\lambda)$ are the wavelength dependent transmission patterns provided by the SLM, which can be calibrated by analyzing its electrically-controlled birefringence. Since the detector array is spatially pixellated, the $(m,n)^{th}$ detector measurement is given by:

where $T_{m,n,j}$ represents the discretized transmission patterns. Fig. 3 is a photograph of the experimental setup. The setup includes a 60 mm objective lens (Jenoptik), a 75 mm achromatic relay lens (Edmound Optics), a polarization beam splitter (Newport), an achromatic quarter wave plate (Newport), a liquid crystal based SLM (Pluto, Holoeye), two 75 mm imaging lenses (Pentax), and a monochrome CCD camera (Pike F-421, AVT) with 2048$\times$2048 pixels that are 7.4 $\mu$m square. A 450-680nm band pass filter (Badder) is mounted on the objective lens to block unwanted ultra-violet and infrared light (the SLM and camera are optimized for visible-range operation). The SLM provides strong modulation and effective multiplexing to fulfill the requirement of compressive sensing. This modulator is based on a reflective Liquid Crystal on Silicon (LCoS) microdisplay technique, which has a 1920$\times$1080-pixel active area with an 8 $\mu$m pixel pitch.

II-B2 System Calibration

The effect of the 8-bit, pseudo-random voltage pattern on the spectral channels can be estimated by recording the system response generated by the SLM. Fig. 4 records the transmission under the monochromatic illumination and homogeneously-applied voltage on the modulator. We average the transmission across the center area of the SLM and calculate its relative transmisson. A potential advantage of this SLM is that we can change the code easily to get different modulations.This also opens the door of adaptive compressive sensing . However, after plenty experiments, we found the contribution of adaptive sensing is very limited in the CS inversion, partially due to the mismatch between the designed code and the calibration. We further found random coding performing well in our camera.

We assume that the system operator provides one-to-one mapping (1:1 magnification of the SLM pixels onto the detector) between the micro-display and the detector array. Theoretically, the system operator $T(x,y,\lambda)$ can be estimated by using SLM’s calibration data and the applied voltage on the SLM. However, one might account the error in the real system projection. For example, optical aberrations and sub-pixel alignment discrepancies between the detector and SLM might break the ideal mapping and image some of the SLM pixels onto several detector pixels. Therefore, part of the transmission code might deviate from the ideal value, which results in an inaccurate system operator. Since the system operator dominates the quality of the object estimation in this inverse problem, a careful calibration of the $T(x,y,\lambda)$ is required.

A better representation of the system response can be acquired by recording the transmission pattern illuminated by each spectral band. The revised system operator $T(x,y,\lambda)$ includes all the possible coding patterns generated by the SLM, which can contribute to the data reconstruction. Here we combine the tungsten light source (LSH-T250, Horiba) with a monochromator (iHR320, Horiba) to quantize the spectral dimension into bands of finite width. Each band has a 7.5-8 nm full width at half maximum (FWHM). The scene’s spectral irradiance is recorded by a fiber optics spectrometer (USB2000, Ocean Optics); these values are taken as ground truth in the experiments presented later in the paper. The number of spectral channels and their central wavelengths are determined by the grating period (1800 periods/mm) and the optical path length of the monochromator. To better represent the continuous light source, two adjacent spectral bands are separated by the monochromator’s FWHM. At this calibration resolution, the system’s spectrally-sensitive bandwidth (450-680nm) has been discretized into $30$ spectral channels. Importantly, the spectral resolution of the reconstructed data is determined by the monochromator’s resolving power during calibration; smaller FWHM values result in larger numbers of calibrated and reconstructed spectral channels. Compared to the spectral imagers that are reliant upon dispersive elements to spatially shear the physical code (in CASSI), this SLM based spectral imager can easily improve the spectral resolution without revising the main camera’s optical design.

II-C Shared Mathematical Model of CASSI and SLM-CASSI

The pixel-wise modulation of CASSI results in feasible patch-based reconstruction algorithms.

Since the SLM-CASSI sensing paradigm reconstructs $N_{\lambda}$ spectral images from a single measurement ${\bf M}$, the compression ratio of the CASSI and SLM-CASSI systems discussed above is $N_{\lambda}$. The compressive sampling process in both spectral imaging systems can be represented by the same matrix described in the following section.

III Reconstruct Hyperspectral Images with Blind Compressive Sensing

In this section, we develop a new blind Bayesian CS model to reconstruct hyperspectral images from CASSI or SLM-CASSI measurements. The proposed model is a generalized framework which can also be used in denoising and inpainting problems, among others .

where $\boldsymbol{\epsilon}_{n}$ represents the additive Gaussian noise, and we model it as ${\boldsymbol{\epsilon}}_{n}\sim{\cal N}(0,\alpha_{0}^{-1}{\bf I}_{P})$, with $\alpha_{0}$ denoting the noise precision.The extension of our model to non-uniform denoising problem is straightforward, i.e., by imposing spatially-varying noise models for different patches. We place a diffuse gamma prior on $\alpha_{0}$

where $(c_{0},d_{0})$ are hyperparameters. Taking account of the dictionary learning model, (7) can be written as:

where ${\boldsymbol{s}_{n}}$ is a vector of coefficients describing the decomposition of the signal $\boldsymbol{x}_{n}$ in terms of dictionary atoms. Given the measurement set ${\bf Y}$ and the forward matrices $\{{\boldsymbol{\Psi}}_{n}\}_{n=1}^{N}$, we aim to jointly learn ${\bf D},\{{\boldsymbol{s}}_{n}\}_{n=1}^{N}$, and the noise precision parameter $\alpha_{0}$ to recover ${\bf X}$. One key difference of our model compared to other blind CS work is that each patch has a unique $\boldsymbol{\Psi}_{n}$, which is inspired from our cameras since the mask is generated randomly and $\boldsymbol{\Psi}_{n}$ also takes account of the system calibration.

We model each dictionary atom as a draw from a Gaussian distribution,

The coefficients $s_{k,n}$ are assumed drawn from the marginalized distribution

where ${\rm InvGa}(\cdot)$ denotes the inverse-gamma distribution. The parameter $\tau_{n}>0$ is a “global” scaling for all coefficients of the $n^{th}$ patch, and $\Phi_{k,n}$ is a “local” weight for the $k$th coefficient of that patch. We place a gamma prior ${\rm Ga}(g_{0},h_{0})$ on $\Phi_{k,n}$, where one may set the hyperparameters $(g_{0},h_{0})$ to ensure that most $\Phi_{k,n}$ are small. This encourages most $s_{k,n}$ to be small. By maximizing the log posterior to obtain a point estimate for the model parameters, one observes that the log of the prior in (III-A) corresponds to adaptive Lasso regularization .

Equivalently to (III-A), the model for $s_{k,n}$ may be represented in the hierarchical form

where latent variables $\{\alpha_{k,n}\}$ are included in the generative model, instead of marginalizing them out as in (III-A). A vague/diffuse gamma prior is placed on the scaling parameters $\tau_{n}$. Despite introducing the latent variables $\{\alpha_{k,n}\}$, the form in (15) is convenient for computation, and with an appropriate choice of $(g_{0},h_{0})$, most $\{\alpha_{k,n}\}$ are encouraged to be large. The large $\alpha_{k,n}$ corresponds to small $s_{k,n}$; this model imposes that most $s_{k,n}$ are small, $i.e.$ it imposes compressibility. Note that (14) is different from the model used in , where a single $\tau$ is used for all patches; here, we impose different compressibility for each patch by inferring a unique $\tau_{n}$, thus providing flexibility.

In order to automatically infer the number of necessary dictionary atoms, we can replace (11) with

The multiplicative gamma prior (MGP) used above is developed to stochastically increase the precision $\eta_{k}$ as $k$ increases. During inference, we observe that as $k$ increases, $\nu_{k}$ tends to zero. This results in an approximate ordering of $\nu_{k}$ by weight. Based on this weight, we can infer the importance of the columns of ${\bf D}$. The other use of the $\nu_{k}$’s is to avoid over-fitting during the learning procedure. We can update a fraction of ${\bf D}$ by the weight of $\nu_{k}$ and discard the atoms with small weights to reduce the computational cost (refer to the Appendix for the inferred $\nu_{k}$ and learned dictionary ${\bf D}$).

III-B The Statistical Model

The full statistical model of the proposed blind CS approach is:

where broad priors are placed on the hyperparameters $(a_{0},b_{0},c_{0},d_{0},e_{0},f_{0},g_{0},h_{0})$, i.e., $a_{0}=\dots=h_{0}=10^{-6}$. Note in (23), the coefficients are also scaled to the noise precision $\alpha_{0}$.

Let ${\bf\Theta}=\{{{\bf D}},{\bf S},{\bf\Lambda},{\boldsymbol{\boldsymbol{\epsilon}}}\}$ denote the parameters to be inferred. The log of the joint posterior may be expressed as:

III-C Related Models

the parameters $\alpha_{0}$ and $\gamma_{s}$ are typically set by hand (e.g., via cross-validation). One advantage of the Bayesian framework is that we infer posterior distributions for $\alpha_{0}$ and $\gamma_{s}$ (in our model this is ($\tau_{n},\alpha_{k,n}$), along with similar posterior estimates for all model parameters without cross-validation. We also note that if we replace (30) by the spike-slab prior as used in , the model will reduce to BPFA. As opposed to this spike-slab prior, which imposes sparsity directly, our shrinkage prior imposes compressibility, and it is more appropriate to the high dimensional hyperspectral image patches. The global-local shrinkage prior used here can extract more useful information from the limited dictionary atomsWe did experiments of denoising and inpainting with benchmark color images and compared with K-SVD and BPFA . The results are shown in the appendix. The proposed algorithm constantly performs better than the above two methods..

Other forms of shrinkage priors like Gaussian scale model and Laplace scale mode can be found in . Aiming to better mimic the marginal behavior of discrete mixture priors, the global-local shrinkage priors have been developed to offer sufficient flexibility in high-dimensional settings, which inspires our model.

III-D Inference

Due to local conjugacy, we can write the conditional posterior distribution for all parameters of our model in closed form, making the following Markov Chain Monte Carlo (MCMC) inference based on Gibbs sampling a straightforward procedure.

where ${\rm IG}(\cdot)$ denotes the inverse-Gaussian distribution.

where ${\rm GIG}(x;a,b,p)$ is the generalized inverse Gaussian distribution

and $K_{p}(\theta)$ is the modified Bessel function of the second kind

where $\|{\boldsymbol{\Psi}}_{n}\|_{0}$ denotes the number of nonzero entries in ${\boldsymbol{\Psi}}_{n}$, and “-” refers to the conditioning parameters of the distributions.

III-E RGB Images as Side Information

While reconstructing the hyperspectral images is challenging (the compression ratio is $N_{\lambda}:1$ when a single measurement is available), an RGB image of the same scene can be obtained easily by an off-the-shelf color camera in the unused path of the system (i.e. directly reflecting off the polarizing beamsplitter) as shown in Fig. 2. We here consider using the RGB image as side information to aid the reconstruction.

In the case of an additional side RGB camera, the measurement is a joint dataset composed of the RGB image and the CASSI measurement; the compression ratio can be considered as $(N_{\lambda}+3):4$. The new measurement model is:

where $\{{\bf Y}^{(\rm rgb)},{\bf X}^{(\rm rgb)}\}$ are the patch format RGB image. Considering each patch,

Another way to use the RGB image is treating each R, G and B channel as a superposition of the hyperspectral images with different weights corresponding to the quantum efficiency of the R, G, and B sensors. However, these quantum efficiencies may be different for each camera. Here, we aim to propose a general and robust framework to collaborate RGB images with CASSI measurements. Therefore, the formulation in (45) is adopted.

In our experimental setup, a grayscale image may also be used as side information. As the RGB camera is inexpensive and carrying richer spectral information, we only consider the RGB case in the following experiments. In the results presented below, we consider that RGB image is measured separately by an additional camera. We are now modifying our SLM-CASSI system to capture the RGB image and the compressed hyperspectral image simultaneously as illustrated in Fig. 2.

IV Experimental Results

We employ a spatial patch size $n_{x}=n_{y}=8$. The compression ratio, $N_{\lambda}$, depends on the dataset. The proposed model can learn the dictionary atom number $K$ from the data. During experimentation, we have found that setting $K=64$ usually provides good results. We have verified that further increasing $K$ does not significantly change the outcome of our methods. All code used in the experiments was written in Matlab and executed on a 3.3GHz desktop with 16GB RAM.

We evaluate the proposed algorithm on synthetic data presented in Section IV-A and the real data captured both by the CASSI and the proposed SLM-CASSI camera Section IV-B. We denote our method as “shrinkage” in the experiments (in figures and ta bles) as the shrinkage prior is used in our model. To evaluate the algorithm’s performance, we use the PSNR of the reconstructed images at different wavelengths and the correlation between the reconstructed spectrum and the true (reference) spectrum.

The computational time of our model is similar to the BPFA model used in and linearized Bregman. TwIST provides faster, but generally worse results than our algorithm in terms of the metrics mentioned above. Quantitatively, linearized Bregman and TwIST require about 20 minutes and 14 minutes, respectively, to reconstruct the bird data (size $768\times 1024\times 24$).

We use hyperspectral images encoded with a random binary (Bernoulli(0.5)) mask to simulate the CASSI measurements. The RGB images are available and aligned well with the hyperspectral images.

We first consider the hyperspectral images of natural scenes used in http://personalpages.manchester.ac.uk/staff/d.h.foster/Hyperspectral_images_of_natural_scenes_04.html. There are $N_{\lambda}=33$ channels (400-720nm with a 10nm interval) and we resize each image to $N_{x}=N_{y}=512$. The PSNR of each reconstructed spectral channel with mean values and standard deviations (across wavelength channels) are presented for each algorithm in Table I. It can be seen that the proposed shrinkage method outperforms the others, especially once the RGB images are used as side information. Fig. 5 plots the reconstructed spectrum of selected blocks (the average spectrum of pixels inside the block is used) in the scene compared to the ground truth. It can be seen that: $i$) though TwIST provides the lowest PSNR, the reconstructed spectrum is usually correct; $ii$) the spectrum reconstructed by the proposed algorithm is improved significantly when side information is provided; $iii$) the linearized Bregman presents the worst spectrum, mainly due to the DCT used in the spectral domain for inversion.

IV-A2 Bird data

Next we consider the 24-channel bird data (Fig. 6) measured by a hyperspectral camera in . The RGB image is also available. Fig. 7 shows the reconstructed images using the proposed algorithm with and without side information. The left part of Fig. 8 plots the reconstruction PSNR of each channel using different algorithms. Again, our shrinkage blind CS method coupled with the RGB image provides the best result. The right part of Fig. 8 shows the reconstructed images of five selected channels using the different algorithms. It can be seen that the TwIST results are characterized by lost details (over smoothed) of the birds; our model’s results without using the RGB images appear noisy. Fig. 9 depicts the reconstructed spectra of different birds with different algorithms. It can be seen that both TwIST and our algorithm with the RGB image provide very good matches to the ground truth, while the proposed model without RGB images and linearized Bregman do not represent the spectrum well.

IV-B Real Data

In this section, we apply the proposed algorithm to real data captured by our cameras (both the original CASSI camera and the new SLM-CASSI camera).

We first demonstrate our algorithm on data taken by the original CASSI system http://www.disp.duke.edu/projects/CASSI/experimentaldata/index.ptml. In these experiments, the reconstructions have 33 spectral channels (marked on the reconstruction in Fig. 10) and $N_{x}=N_{y}=256$. There are 4 objects in the scene, a red apple, a yellow banana, a blue stapler and a green pineapple. Fig. 10 shows the reconstruction of the proposed algorithm without the RGB image. Since the RGB image is not well-aligned with the CASSI measurement, we don’t show the result with side information and we found the reconstruction with the RGB image (not aligned) looks similar to this one due to the simple scene used in this experiment. Fig. 11 compares selected reconstructed images inverted by different algorithms. It can be seen that the proposed algorithm provides more detail of the scene (notice the clear apple stem). Fig. 12 plots the reconstructed spectra of the four objects. TwIST has yielded accurate spectra; our algorithm (without side information) provides similar results. The linearized Bregman algorithm does not reconstruct this real data well; we don’t not show its results in the following experiments.

IV-B2 Bird data

We again consider the bird data, now captured by the multiframe-CASSI camera . The RGB image is aligned manually with the CASSI measurement ($N_{x}=703,N_{y}=1021,N_{\lambda}=24$). We plot the spectra in Fig. 13; the reconstructed images are shown in Fig. 14. It can be seen clearly that our proposed method with side information provides the best results with respect to both image clarity and spectral accuracy. Without the use of side information, the proposed model yields clear images but reconstructs the spectra poorly. This verifies the benefit of the side information in real data.

IV-B3 Data with CASSI-SLM

Now we consider the dataset captured by the proposed camera. Since no RGB image is available, we only show the results of our algorithm without side information. Fig. 15 shows the reconstructed spectrum and selected frames for the M&M dataset ($N_{x}=512,N_{y}=784,N_{\lambda}=30$)The 30 wavelengths are 450nm, 458nm, 465nm, 473nm, 481.5nm, 489.5nm, 498nm, 507nm, 516nm, 524.5nm, 532.5nm, 540.5nm, 548.5nm, 556,5nm, 564.5nm, 572.5nm, 580.5nm, 588.5nm, 596nm, 603.5nm, 611nm, 618.5nm, 625.8nm, 633.5nm, 641nm, 648.5nm, 656nm, 663.5nm, 671nm, 678.5nm.. As with the original CASSI experiments, we use a fiber optic spectrometer (USB2000, Ocean Optics) to provide the reference spectra for the targets. It can be seen again our algorithm provides better results than TwIST, both for images and spectrum.

As a last example, we show the reconstructed images with our method for the berry data ($N_{x}=N_{y}=2048,N_{\lambda}=40$) in Figure 16. Notice that the leaf reconstructs are prevalent in the green channels (520$\sim$590nm), while the berries (red) appear in the red portion of the spectrum (610$\sim$680nm).

V Conclusions

We have developed and tested a new Bayesian dictionary learning model for blind CS. Specifically, we have demonstrated high-quality inversion of compressed hyperspectral images captured by real cameras. Via the global-local shrinkage priors, our algorithm imposes compressibility on the dictionary coefficients to extract more information from the dictionary atoms. The reconstruction quality is improved significantly by integrating the compressed measurements with RGB images as side information under the blind CS framework. We also developed a new compressive hyperspectral imaging camera that uses an SLM to perform the spectral coding. Experimental results demonstrate the feasibility of the new camera, the superior performance of the algorithm, and the benefit of side information.

The original CASSI camera utilizes a separate disperser and encoder; the SLM-CASSI is capable of modulating both dimensions with one element. The SLM-CASSI also provides the potential of video rate compressive sensing and multi-frames adaptive sensing using the 60 Hz refresh rate of the SLM. SLM-CASSI offers greater flexibility of masking functions and an additional RGB camera for side information; however, it intertwines spectral and spatial multiplexing capabilities. The coded aperture/disperser architecture used in the original CASSI camera is more reliable, has a smaller form factor (i.e. fewer optical parts), and is less expensive, but lacks the coding flexibility and ability to use a separate camera for side information.

In this paper, we have shown that by using the RGB image as side information in our proposed blind CS framework, the reconstruction quality of hyperspectral images can be improved significantly. This framework can be generalized to other applications and algorithms. For instance, Gaussian mixture model based dictionary learning approaches [44, 45, 46, Yang14GMM2, 47] that learn a union of subspaces can also benefit from side information.

References

VI Model and Inference

The posterior density function of model parameters may be represented as:

where ${\rm IG}$ denotes inverse-Gaussian distribution.

where ${\rm GIG}(x:a,b,p)$ is the generalized inverse Gaussian distribution:

and $K_{p}(\theta)$ is the modified Bessel function of the second kind

where $\|{\boldsymbol{\Psi}}_{n}\|_{0}$ denotes the number of nonzero entries in ${\boldsymbol{\Psi}}_{n}$.

VI-B Variational Bayesian Inference

We adopt a mean-field variational Bayesian (VB) inference in lieu of its improved runtime. A VB approach attempts to approximate the posterior distribution by a simpler distribution, $p({\bf\Theta}|{\boldsymbol{Y}})\approx q({\bf\Theta})$, where ${\bf Y}$ is the observed data matrix and ${\bf\Theta}$ denotes the set of independent latent variables in the model . VB assumes a complete factorization across latent variables, $q({\bf\Theta})=\prod_{i}q_{i}({\bf\Theta}_{i})$. We define ${\bf\Theta}=\{{\boldsymbol{D}},{\boldsymbol{S}},{\boldsymbol{\Lambda}},{\boldsymbol{\boldsymbol{\epsilon}}}\}$ for the purpose of this work.

Solving for the optimal distribution $q^{\star}({\bf\Theta})$ that minimizes the distance between $p$ and $q$ effectively estimates the conditional posterior distribution $p({\bf\Theta}|{\bf Y})$. A commonly-used distance metric between the two distributions functions is the Kullback-Leibler (KL) divergence . We write the KL-divergence of $p$ from $q$ as follows:

We here observe that $\ln p({\bf Y})$ is fixed with respect to the variations in $q({\bf\Theta})$. Therefore, maximizing the Evidence Lower Bound (ELBO) $\L(q)$ is equivalent to minimizing the KL-divergence between the two distributions. This minimal distance occurs when

Assuming a complete factorization across the latent variables $q({\bf\Theta})=\prod_{i}q_{i}({\bf\Theta}_{i})$, each parameter in a variational Bayes model is independently updated according to

The update equations of VB are straightforward from the posterior distribution of Gibbs sampling ($\langle\cdot\rangle$ represents the expectation of the random variable inside):

with ${\rm Tr(\cdot)}$ denoting the trace of the matrix inside $(\hskip 0.85358pt)$.

VII More Results

To demonstrate the general applications of our proposed dictionary learning model, we show the denoising and inpainting results from benchmark color images (8-bits) tested in . For inpainting, we show PSNRs of the restored images at various observed data ratios ($20\%$ means $80\%$ pixels are missing) in Table II. Noisy images corrupted with zero-mean Gaussian noise with different standard deviations $\sigma$ are restored with PSNRs shown in Table III. Our proposed algorithm consistently provides the highest PSNR. Inpainting results for corrupted images (Figure 18) are presented in Figure 19 (PSNRs shown in Table II); denoising results for noisy images (corrupted with zero-mean Gaussian noise of various standard deviations, Figure 20) are presented in Figure 21 (PSNRs shown in Table III). Figure 17 shows an example of a dictionary and $\nu_{k}$ learned from a noisy image using the proposed dictionary learning model (Section VI). Importantly, few iterations of our VB inference are required to obtain good results; 20 iterations are used for denoising and 100 iterations are used for inpainting ($\sim$1 second per iteration with an image size $256\times 256\times 3$ on an i5 CPU with non-optimized MATLAB code).