Compressive Hyperspectral Imaging via Approximate Message Passing

I Introduction

A hyperspectral image is a three-dimensional (3D) image cube comprised of a collection of two-dimensional (2D) images (slices), where each 2D image is captured at a specific wavelength. Hyperspectral images allow us to analyze spectral information about each spatial point in a scene, and thus can help us identify different materials that appear in the scene . Therefore, hyperspectral imaging has applications to areas such as medical imaging , remote sensing , geology , and astronomy .

Conventional spectral imagers include whisk broom scanners, push broom scanners , and spectrometers . In whisk broom scanners, a mirror reflects light onto a single detector, so that one pixel of data is collected at a time; in push broom scanners, an image cube is captured with one focal plane array (FPA) measurement per spatial line of the scene; and in spectrometers, a set of optical bandpass filters are tuned in steps in order to scan the scene. The disadvantages of these techniques are that (i) data acquisition takes a long time, because they require scanning a number of zones linearly in proportion to the desired spatial and spectral resolution; and (ii) large amounts of data are acquired and must be stored and transmitted. For example, for a megapixel camera ($10^{6}$ pixels) that captures a few hundred spectral bands ($>100$ spectral channels) at 8 or 16 bits per frame, conventional spectral imagers demand roughly 10 megabytes per raw spectral image, and thus require space on the order of gigabytes for transmission or storage, which exceeds existing streaming capabilities.

To address the limitations of conventional spectral imaging techniques, many spectral imager sampling schemes based on compressive sensing have been proposed . The coded aperture snapshot spectral imager (CASSI) is a popular compressive spectral imager and acquires image data from different wavelengths simultaneously. In CASSI, the voxels of a scene are first coded by an aperture, then dispersed by a dispersive element, and finally detected by a 2D FPA. That is, a 3D image cube is suppressed and measured by a 2D array, and thus CASSI acquires far fewer measurements than those acquired by conventional spectral imagers, which significantly accelerates the imaging process. In particular, for a data cube with spatial resolution of $N\times M$ and $L$ spectral bands, conventional spectral imagers collect $MNL$ measurements. In contrast, CASSI collects measurements on the order of $M(N+L-1)$. Therefore, the acquisition time, storage space, and required bandwidth for transmission in CASSI are reduced. On the other hand, because the measurements from CASSI are highly compressive, reconstructing 3D image cubes from CASSI measurements becomes challenging. Moreover, because of the massive size of 3D image data, it is desirable to develop fast reconstruction algorithms in order to realize real time acquisition and processing.

Fortunately, it is possible to reconstruct the 3D cube from the 2D measurements according to the theory of compressive sensing , because the 2D images from different wavelengths are highly correlated, and the 3D image cube is sparse in an appropriate transform domain, meaning that only a small portion of the transform coefficients have large values. Approximate message passing (AMP) has recently become a popular algorithm that solves compressive sensing problems, owing to its promising performance and efficiency. Therefore, we are motivated to investigate how to apply AMP to the CASSI system.

I-B Related work

Despite the good results attained with the algorithms mentioned above, they all need manual tuning of some parameters, which may be time consuming. In GPSR and TwIST, the optimal regularization parameter could be different in reconstructing different image cubes. In dictionary learning methods, although the parameters can be learned automatically by methods such as Markov Chain Monte Carlo, the learning process is usually time consuming. Moreover, the patch size and the number of dictionary atoms in dictionary learning methods must be chosen carefully.

I-C Contributions

In this paper, we develop a robust and fast reconstruction algorithm for the CASSI system using approximate message passing (AMP) . AMP is an iterative algorithm that can apply image denoising at each iteration. Previously, we proposed a 2D compressive imaging reconstruction algorithm, AMP-Wiener , where an adaptive Wiener filter was applied as the image denoiser within AMP. Our numerical results showed that AMP-Wiener outperformed the prior art in terms of both reconstruction quality and runtime. The current paper extends AMP-Wiener to reconstruct 3D hyperspectral images from the CASSI system, and we call the new approach “AMP-3D-Wiener.” Because the matrix that models the CASSI system is highly sparse, structured, and ill-conditioned, applying AMP to the CASSI system becomes challenging. For example, (i) the noisy image cube that is obtained at each AMP iteration contains non-Gaussian noise; and (ii) AMP encounters divergence problems, i.e., the reconstruction error may increase with more iterations. Although it is favorable to use a high-quality denoiser within AMP, so that the reconstruction error may decrease faster as the number of iteration increases, we have found that in such an ill-conditioned imaging system, applying aggressive denoisers within AMP causes divergence problems. Therefore, besides using standard techniques such as damping to encourage the convergence of AMP, we modify the adaptive Wiener filter and make it robust to the ill-conditioned system model. There are existing denoisers that may outperform the modified adaptive Wiener filter in a single step denoising problem. However, the modified adaptive Wiener filter fits into the AMP framework and allows AMP to improve over successive iterations.

Our approach is applied in nature, and the convergence of AMP-3D-Wiener is tested numerically. We simulate AMP-3D-Wiener on several settings where complementary random coded apertures (see details in Section IV-A) are employed. The numerical results show that AMP-3D-Wiener reconstructs 3D image cubes with less runtime and higher quality than other compressive hyperspectral imaging reconstruction algorithms such as GPSR and TwIST (Figure 3), even when the regularization parameters in GPSR and TwIST have already been tuned. These favorable results provide AMP-3D-Wiener major advantages over GPSR and TwIST. First, when the bottleneck is the time required to run the reconstruction algorithm, AMP-3D-Wiener can provide the same reconstruction quality in 100 seconds that the other algorithms provide in 450 seconds (Figure 3). Second, when the bottleneck is the time required for signal acquisition by CASSI hardware, the improved reconstruction quality could allow to reduce the number of shots taken by CASSI by as much as a factor of $2$ (Figure 8). Finally, the reconstructed image cube can be obtained by running AMP-3D-Wiener only once, because AMP-3D-Wiener does not need to tune any parameters. In contrast, the regularization parameters in GPSR and TwIST need to be tuned carefully, because the optimal values of these parameters may vary for different test image cubes. In order to tune the parameters for each test image cube, we run GPSR and TwIST many times with different parameter values, and then select the ones that provide the best results.

The remainder of the paper is arranged as follows. We review CASSI in Section II, and describe our AMP based compressive hyperspectral imaging reconstruction algorithm in Section III. Numerical results are presented in Section IV, while Section V concludes.

II Coded Aperture Snapshot Spectral Imager (CASSI)

The coded aperture snapshot spectral imager (CASSI) is a compressive spectral imaging system that collects far fewer measurements than traditional spectrometers. In CASSI, (i) the 2D spatial information of a scene is coded by an aperture, (ii) the coded spatial projections are spectrally shifted by a dispersive element, and (iii) the coded and shifted projections are detected by a 2D FPA. That is, in each coordinate of the FPA, the received projection is an integration of the coded and shifted voxels over all spectral bands at the same spatial coordinate. More specifically, let $f_{0}(x,y,\lambda)$ denote the voxel intensity of a scene at spatial coordinate $(x,y)$ and at wavelength $\lambda$, and let $T(x,y)$ denote the coded aperture. The coded density $T(x,y)f_{0}(x,y,\lambda)$ is then spectrally shifted by the dispersive element along one of the spatial dimensions. The energy received by the FPA at coordinate $(x,y)$ is therefore

where $S(\lambda)$ is the dispersion function induced by the prism at wavelength $\lambda$. Suppose we take a scene of spatial dimension $M$ by $N$ and spectral dimension $L$, i.e., the dimension of the image cube is $M\times N\times L$, and the dispersion is along the second spatial dimension $y$, then the number of measurements captured by the FPA will be $M(N+L-1)$. If we approximate the integral in (1) by a discrete summation and vectorize the 3D image cube and the 2D measurements, then we obtain a matrix-vector form of (1),

where ${\bf f_{0}}$ is the vectorized 3D image cube of dimension $n=MNL$, vectors ${\bf g}$ and ${\bf z}$ are the measurements and the additive noise, respectively, and the matrix ${\bf H}$ is an equivalent linear operator that models the integral in (1). In this paper, we assume that the additive noise ${\bf z}$ is independent and identically distributed (i.i.d.) Gaussian. With a single shot of CASSI, the number of measurements is $m=M(N+L-1)$, whereas $K$ shots will yield $m=KM(N+L-1)$ measurements. The matrix H in (2) accounts for the effects of the coded aperture and the dispersive element. A sketch of this matrix is depicted in Figure 1(a) when $K=2$ shots are used. It consists of a set of diagonal patterns that repeat in the horizontal direction, each time with a unit downward shift, as many times as the number of spectral bands. Each diagonal pattern is the coded aperture itself after being column-wise vectorized. Just below, the next set of diagonal patterns are determined by the coded aperture pattern used in the subsequent shot. The matrix H will thus have as many sets of diagonal patterns as FPA measurements. Although ${\bf H}$ is sparse and highly structured, the restricted isometry property still holds, as shown by Arguello and Arce .

II-B Higher order CASSI

Recently, Arguello et al. proposed a higher order model to characterize the CASSI system with greater precision, and improved the quality of the reconstructed 3D image cubes. In the standard CASSI system model, each cubic voxel in the 3D cube contributes to exactly one measurement in the FPA. In the higher order CASSI model, however, each cubic voxel is shifted to an oblique voxel because of the continuous nature of the dispersion, and therefore the oblique voxel contributes to more than one measurement in the FPA. As a result, the matrix ${\bf H}$ in (2) will have multiple diagonals as shown in Figure 1(b), where there are sets of $3$ diagonals for each FPA shot, accounting for the voxel energy impinging into the neighboring FPA pixels. In this case, the number of measurements with $K=1$ shot of CASSI will be $m=M(N+L+1)$, because each diagonal entails the use of $M$ more pixels (we refer readers to for details).

In Section IV, we will provide promising image reconstruction results for this higher order CASSI system. Using the standard CASSI model, our proposed algorithm produces similar advantageous results over other competing algorithms.

III Proposed Algorithm

The goal of our proposed algorithm is to reconstruct the image cube ${\bf f_{0}}$ from its compressive measurements ${\bf g}$, where the matrix ${\bf H}$ is known. In this section, we describe our algorithm in detail. The algorithm employs (i) approximate message passing (AMP) , an iterative algorithm for compressive sensing problems, and (ii) adaptive Wiener filtering, a hyperspectral image denoiser that can be applied within each iteration of AMP.

Below we describe that the linear imaging system model in (2) can be converted to a 3D image denoising problem in scalar channels. Therefore, we begin by defining scalar channels, where the noisy observations ${\bf q}$ of the image cube ${\bf f_{0}}$ obey

and ${\bf v}$ is the additive noise vector. Recovering ${\bf f_{0}}$ from ${\bf q}$ is known as a 3D image denoising problem.

III-B Approximate message passing

Algorithm framework: AMP has recently become a popular algorithm for solving signal reconstruction problems in linear systems as defined in (2). The AMP algorithm proceeds iteratively according to

where ${\bf H}^{T}$ is the transpose of ${\bf H}$, $R=m/n$ represents the measurement rate, $\eta_{t}(\cdot)$ is a denoising function at the $t$-th iteration, $\eta_{t}^{\prime}({\bf s})=\frac{\partial}{\partial{\bf s}}\eta_{t}({\bf s})$, and $\langle{\bf u}\rangle=\frac{1}{n}\sum_{i=1}^{n}u_{i}$ for some vector ${\bf u}=(u_{1},u_{2},\ldots,u_{n})$. We will explain in Section III-E how ${\bf f}^{t}$ and ${\bf r}^{t}$ are initialized. The last term in (5) is called the “Onsager reaction term” in statistical physics. This Onsager reaction term helps improve the phase transition (trade-off between the measurement rate and signal sparsity) of the reconstruction process over existing iterative thresholding algorithms . In the $t$-th iteration, we obtain the estimated image cube ${\bf f}^{t}$ and the residual ${\bf r}^{t}$. We highlight that the vector ${\bf H}^{T}{\bf r}^{t}+{\bf f}^{t}$ in (4) can be regarded as a noise-corrupted version of ${\bf f_{0}}$ in the $t$-th iteration with noise variance $\sigma_{t}^{2}$, and therefore $\eta_{t}(\cdot)$ is a 3D image denoising function that is performed on a scalar channel as in (3). Let us denote the equivalent scalar channel at iteration $t$ by

where the noise level $\sigma^{2}_{t}$ is estimated by ,

and $r^{t}_{i}$ denotes the $i$-th component of the vector ${\bf r}^{t}$ in (5).

Moreover, if the matrix is i.i.d. random, then the noise in the scalar channel (3) can be viewed as asymptotically i.i.d. Gaussian .

III-C Damping

We have discussed in Section III-B that many mathematical properties of AMP hold for the setting where the measurement matrix is i.i.d. Gaussian. When the measurement matrix is not i.i.d. Gaussian, such as the highly structured matrix ${\bf H}$ defined in (2), AMP may encounter divergence issues. A standard technique called “damping” is frequently employed to solve for the divergence problems of AMP, because it only increases the runtime modestly.

Specifically, damping is an extra step within AMP iterations. In (4), instead of updating the value of ${\bf f}^{t+1}$ by the output of the denoiser $\eta_{t}({\bf H}^{T}{\bf r}^{t}+{\bf f}^{t})$, we assign a weighted average of $\eta_{t}({\bf H}^{T}{\bf r}^{t}+{\bf f}^{t})$ and ${\bf f}^{t}$ to ${\bf f}^{t+1}$ as follows,

for some constant $0<\alpha\leq 1$. Similarly, after obtaining ${\bf r}^{t}$ in (5), we add an extra damping step that updates the value of ${\bf r}^{t}$ to be $\alpha\cdot{\bf r}^{t}+(1-\alpha)\cdot{\bf r}^{t-1}$, where the value of $\alpha$ is the same as that in (8).

AMP has been proved to converge with sufficient damping, under the assumption that the prior of ${\bf f_{0}}$ is i.i.d. Gaussian with fixed means and variances throughout all iterations, and the amount of damping depends on the condition number of the matrix ${\bf H}$. Note that other AMP variants have also been proposed in order to encourage convergence for a broader class of measurement matrices.

in our modified algorithm AMP-3D-Wiener, we propose a simpler version of adaptive Wiener filter as described in Section III-D to stabilize the estimation of the prior distribution of ${\bf f_{0}}$. Although we do not have justifications for convergence of AMP-3D-Wiener at this point, we find in our simulations that AMP-3D-Wiener converges for all tested hyperspectral image cubes with moderate amount of damping.

III-D Adaptive Wiener filter

We are now ready to describe our 3D image denoiser, which is the function $\eta_{t}(\cdot)$ in the first step of AMP iterations in (4).

Sparsifying transform: Recall that in 2D image denoising problems, a 2D wavelet transform is often performed, and some shrinkage function is applied to the wavelet coefficients in order to suppress noise . The wavelet transform based image denoising method is effective, because natural images are usually sparse in the wavelet transform domain, i.e., there are only a few large wavelet coefficients and the rest of the coefficients are small. Therefore, large wavelet coefficients are likely to contain information about the image, whereas small coefficients are usually comprised mostly of noise, and so it is effective to denoise by shrinking the small coefficients toward zero and suppressing the large coefficients according to the noise variance. Similarly, in hyperspectral image denoising, we want to find a sparsifying transform such that hyperspectral images have only a few large coefficients in this transform domain. Inspired by Arguello and Arce , we apply a wavelet transform to each of the 2D images in a 3D cube, and then apply a discrete cosine transform (DCT) along the spectral dimension, because the 2D slices from different wavelengths are highly correlated. That is, the sparsifying transform ${\bf\Psi}$ can be expressed as a Kronecker product of a DCT transform ${\bf\Phi}$ and a 2D wavelet transform ${\bf W}$, i.e., ${\bf\Psi=\Phi\otimes W}$, and it can be shown that ${\bf\Psi}$ is an orthonormal transform. Let ${\bf\theta}_{\bf q}^{t}$ denote the coefficients of ${\bf q}^{t}$ in this transform domain, i.e., ${\bf\theta}_{{\bf q}}^{t}={\bf\Psi q}^{t}$. Our 3D image denoising procedure will be applied to the coefficients ${\bf\theta}_{\bf q}^{t}$. Besides 2D wavelet transform and 1D DCT, it is also possible to sparsify 3D image cubes by dictionary learning or Gaussian mixture models . Moreover, using an endmember mixing matrix is an alternative to DCT for characterizing the spectral correlation of 3D image cubes. In this work, we focus on a 2D wavelet transform and 1D DCT as the sparsifying transform, because it is an efficient transform that does not depend on any particular types of image cubes, and an orthonormal transform that is suitable for the AMP framework.

Parameter estimation in the Wiener filter: In our previous work on compressive imaging reconstruction problems for 2D images, one of the image denoisers we employed was an adaptive Wiener filter in the wavelet domain, where the variance of each wavelet coefficient was estimated from its neighboring coefficients within a $5\times 5$ window, i.e., the variance was estimated locally.

As an initial attempt, we applied the previously proposed AMP-Wiener to the reconstruction problem in the CASSI system defined in (2). More specifically, the previously proposed adaptive Wiener filter is applied to the noisy coefficients $\theta_{\bf q}^{t}$. Unfortunately, AMP-Wiener encounters divergence issues for the CASSI system even with significant damping such as $\alpha=0.01$ in (8). AMP-Wiener diverges, because it is designed for the setting where the measurement matrix is i.i.d. Gaussian, whereas the measurement matrix ${\bf H}$ defined in (2) is highly structured and not i.i.d., and we found in our numerical experiments that the scalar channel noise ${\bf v}^{t}$ in (6) is not i.i.d. Gaussian. On the other hand, because the Wiener filter allows to conveniently calculate the Onsager term in (5), we are motivated to keep the Wiener filter strategy, although the scalar channel (6) does not contain i.i.d. Gaussian noise. Seeing that estimating the coefficient variance from its neighboring coefficients (a $3\times 3$ or $5\times 5$ neighboring window) does not produce reasonable reconstruction for the CASSI system, we modify the local variance estimation to a global estimation within each wavelet subband. The coefficients $\widehat{\bf\theta}_{\bf f}^{t}$ of the estimated (denoised) image cube ${\bf f}^{t}$ are obtained by Wiener filtering, which can be interpreted as the conditional expectation of ${\bf\theta}_{\bf f}$ given ${\bf\theta}_{\bf q}^{t}$ under the assumption of Gaussian prior and Gaussian noise,

where ${\bf\theta}_{{\bf q},i}^{t}$ is the $i$-th element of ${\bf\theta}_{\bf q}^{t}$, and $\widehat{\mu}_{i,t}$ and $\widehat{\nu}_{i,t}^{2}$ are the empirical mean and variance of ${\bf\theta}_{{\bf q},i}^{t}$ within an appropriate wavelet subband, respectively. Taking the maximum between 0 and $(\widehat{\nu}_{i,t}^{2}-\widehat{\sigma}_{t}^{2})$ ensures that if the empirical variance $\widehat{\nu}_{i,t}^{2}$ of the noisy coefficients is smaller than the noise variance $\widehat{\sigma}_{t}^{2}$, then the corresponding noisy coefficients are set to 0. After obtaining the denoised coefficients $\widehat{\bf\theta}_{\bf f}^{t}$, the estimated image cube in the $t$-th iteration satisfies ${\bf f}^{t}={\bf\Psi}^{-1}\widehat{\bf\theta}_{\bf f}^{t}={\bf\Psi}^{T}\widehat{\bf\theta}_{\bf f}^{t}$. Therefore, the adaptive Wiener filter as a denoiser function $\eta_{t}(\cdot)$ can be written as

where 0 is a zero matrix, $\widehat{\bf V}_{t}$ is a diagonal matrix with $\widehat{\nu}_{i,t}^{2}$ on its diagonal, ${\bf I}$ is the identify matrix, $\widehat{\boldsymbol{\mu}}_{t}$ is a vector that contains $\widehat{\mu}_{i,t}$, and $\max\{\cdot,\cdot\}$ is operating entry-wise.

We apply this modified adaptive Wiener filter within AMP, and call the algorithm “AMP-3D-Wiener.” We will show in Section IV that only a moderate amount of damping is needed for AMP-3D-Wiener to converge.

III-E Derivative of adaptive Wiener filter

The adaptive Wiener filter described in Section III-D is applied in (4) as the 3D image denoising function $\eta_{t}(\cdot)$. The following step in (5) requires $\eta^{\prime}_{t}(\cdot)$, i.e., the derivative of $\eta_{t}(\cdot)$. We now show how to obtain $\eta^{\prime}_{t}(\cdot)$. It has been discussed that when the sparsifying transform is orthonormal, the derivative calculated in the transform domain is equivalent to the derivative in the image domain. According to (9), the derivative of the Wiener filter in the transform domain with respect to $\widehat{\bf\theta}_{{\bf q},i}^{t}$ is $\max\{0,\widehat{\nu}_{i,t}^{2}-\widehat{\sigma}_{t}^{2}\}/\widehat{\nu}_{i,t}^{2}$. Because the sparsifying transform ${\bf\Psi}$ is orthonormal, the Onsager term in (5) can be calculated efficiently as

where $\mathcal{I}$ is the index set of all image cube elements, and the cardinality of $\mathcal{I}$ is $n=MNL$.

We focus on image denoising in an orthonormal transform domain and apply Wiener filtering to suppress noise, because it is convenient to obtain the Onsager correction term in (5). On the other hand, other denoisers that are not wavelet-DCT based can also be applied within the AMP framework. Metzler et al. , for example, proposed to utilize a block matching and 3D filtering denoising scheme (BM3D) within AMP for 2D compressive imaging reconstruction, and run Monte Carlo to approximate the Onsager correction term. However, the Monte Carlo technique is accurate only when the scalar channel (6) is Gaussian. In the CASSI system model (2), BM4D may be an option for the 3D image denoising procedure. However, because the matrix ${\bf H}$ is ill-conditioned, the scalar channel (6) that is produced by AMP iterations (4,5) is not Gaussian, and thus the Monte Carlo technique fails to approximate the Onsager correction term.

Having completed the description of AMP-3D-Wiener, we summarize AMP-3D-Wiener in Algorithm 1, where $\widehat{\bf f}_{\text{AMP}}$ denotes the image cube reconstructed by AMP-3D-Wiener. Note that in the first iteration of Algorithm 1, initialization of ${\bf q}^{0}$ and $\widehat{\sigma}^{2}_{0}$ may not be necessary, because ${\bf r}^{0}$ is an all-zero vector, and the Onsager term is 0 at iteration 1.

IV Numerical Results

In this section, we provide numerical results where we compare the reconstruction quality and runtime of AMP-3D-Wiener, gradient projection for sparse reconstruction (GPSR) , and two-step iterative shrinkage/thresholding (TwIST) . In all experiments, we use the same coded aperture pattern for AMP-3D-Wiener, GPSR, and TwIST. In order to quantify the reconstruction quality of each algorithm, the peak signal to noise ratio (PSNR) of each 2D slice in the reconstructed cubes is measured. The PSNR is defined as the ratio between the maximum squared value of the ground truth image cube ${\bf f_{0}}$ and the mean square error of the estimation $\widehat{\bf f}$, i.e.,

where $f_{(x,y,\lambda)}$ denotes the element in the cube ${\bf f}$ at spatial coordinate $(x,y)$ and spectral coordinate $\lambda$.

In AMP, the damping parameter $\alpha$ is set to be 0.2. Recall that increasing the amount of damping helps prevent the divergence of AMP-3D-Wiener, and that the divergence issue can be identified by evaluating the values of $\widehat{\sigma}_{t}^{2}$ from (7). We select 0.2 as the damping parameter value, because 0.2 is the maximum damping value such that AMP-3D-Wiener converges in all the image cubes we test. The divergence issues of AMP-3D-Wiener can be detected by evaluating the value of $\widehat{\sigma}_{t}^{2}$ obtained by (7) as a function of iteration number $t$. Recall that $\widehat{\sigma}_{t}^{2}$ estimates the amount of noise in the noisy image cube ${\bf q}^{t}$ at iteration $t$. If AMP-3D-Wiener converges, then we expect the value of $\widehat{\sigma}_{t}^{2}$ to decrease as $t$ increases. Otherwise, we know that AMP-3D-Wiener diverges. The choice of damping mainly depends on the structure of the imaging model in (2) but not on the characteristics of the image cubes, and thus the value of the damping parameter $\alpha$ need not be tuned in our experiments.

To reconstruct the image cube ${\bf f_{0}}$, GPSR and TwIST minimize objective functions of the form

where $\phi(\cdot)$ is a regularization function that characterizes the structure of the image cube ${\bf f_{0}}$, and $\beta$ is a regularization parameter that balances the weights of the two terms in the objective function. In GPSR, $\phi({\bf f})=\|{\bf\Psi f}\|_{1}$; in TwIST, the total variation regularizer is employed,

The first set of simulations is performed for the scene shown in Figure 2. This data cube was acquired using a wide-band Xenon lamp as the illumination source, modulated by a visible monochromator spanning the spectral range between $448$ nm and $664$ nm, and each spectral band has $9$ nm width. The image intensity was captured using a grayscale CCD camera, with pixel size $9.9$ $\mu$m, and 8 bits of intensity levels. The resulting test data cube has $M\times N=256\times 256$ pixels of spatial resolution and $L=24$ spectral bands.

Setting 1: The measurements ${\bf g}$ are captured with $K=2$ shots. The coded aperture in the first shot is generated randomly with 50% of the aperture being opaque, and the coded aperture in the second shot is the complement of the aperture in the first shot. The measurement rate with two shots is $m/n=KM(N+L+1)/(MNL)\approx 0.09$. Moreover, we add Gaussian noise with zero mean to the measurements. The signal to noise ratio (SNR) is defined as $10\log_{10}(\mu_{g}/\sigma_{\text{noise}})$ , where $\mu_{g}$ is the mean value of the measurements ${\bf Hf_{0}}$ and $\sigma_{\text{noise}}$ is the standard deviation of the additive noise ${\bf z}$. In Setting 1, we add measurement noise such that the SNR is 20 dB.

We note in passing that the complementary random coded apertures are binary, and can be implemented through photomask technology or emulated by a digital micromirror device (DMD). Therefore, the complementary random coded apertures are feasible in practice . Moreover, the complementary random coded apertures ensure that in the matrix ${\bf H}$ in (2), the norm of each column is similar, which is suitable for the AMP framework. However, it is a limitation of the current AMP-3D-Wiener that the complementary random coded apertures must be employed, otherwise, AMP-3D-Wiener may diverge.

Let us now evaluate the numerical results for Setting 1. Figure 3 compares the reconstruction quality of AMP-3D-Wiener, GPSR, and TwIST within a certain amount of runtime. Runtime is measured on a Dell OPTIPLEX 9010 running an Intel(R) CoreTM i7-860 with 16GB RAM, and the environment is Matlab R2013a. In Figure 3, the horizontal axis represents runtime in seconds, and the vertical axis is the averaged PSNR over the 24 spectral bands. Although the PSNR of AMP-3D-Wiener oscillates during the first few iterations, which may be because the matrix ${\bf H}$ is ill-conditioned, it becomes stable after 50 seconds and reaches a higher level when compared to the PSNRs of GPSR and TwIST at 50 seconds. After 450 seconds, the average PSNR of the cube reconstructed by AMP-3D-Wiener (solid curve with triangle markers) is 26.16 dB, while the average PSNRs of GPSR (dash curve with circle markers) and TwIST (dash-dotted curve with cross markers) are 23.46 dB and 25.10 dB, respectively. Note that in 450 seconds, TwIST runs roughly 200 iterations, while AMP-3D-Wiener and GPSR run 400 iterations.

Figure 4 complements Figure 3 by illustrating the PSNR of each 2D slice in the reconstructed cube separately. It is shown that the cube reconstructed by AMP-3D-Wiener has $2-4$ dB higher PSNR than the cubes reconstructed by GPSR and $0.4-3$ dB higher than those of TwIST for all 24 slices.

In Figure 5, we plot the 2D slices at wavelengths $488$ nm, $533$ nm, and $578$ nm in the actual image cubes reconstructed by AMP-3D-Wiener, GPSR, and TwIST. The images in these four rows are slices from the ground truth image cube ${\bf f_{0}}$, the cubes reconstructed by AMP-3D-Wiener, GPSR, and TwIST, respectively. The images in columns $1-3$ show the upper-left part of the scene, whereas images in columns $4-6$ show the upper-right part of the scene. All images are of size $128\times 128$. It is clear from Figure 5 that the 2D slices reconstructed by AMP-3D-Wiener have better visual quality; the slices reconstructed by GPSR have blurry edges, and the slices reconstructed by TwIST lack details, because the total variation regularization tends to constrain the images to be piecewise constant.

Furthermore, a spectral signature plot analyzes how the pixel values change along the spectral dimension at a fixed spatial location, and we present such spectral signature plots for the image cubes reconstructed by AMP-3D-Wiener, GPSR, and TwIST in Figure 6. Three spatial locations are selected as shown in Figure 6(a), and the spectral signature plots for locations B, C, and D are shown in Figures 6(b)–6(d), respectively. It can be seen that the spectral signatures of the cube reconstructed by AMP-3D-Wiener closely resemble those of the ground truth image cube (dotted curve with square markers), whereas there are discrepancies between the spectral signatures of the cube reconstructed by GPSR or TwIST and those of the ground truth cube.

According to the runtime experiment from Setting 1, we run AMP-3D-Wiener with 400 iterations, GPSR with 400 iterations, and TwIST with 200 iterations for the rest of the simulations, so that all algorithms complete within a similar amount of time.

Setting 2: In this experiment, we add measurement noise such that the SNR varies from 15 dB to 40 dB, which is the same setting as in Arguello and Arce , and the result is shown in Figure 7. Again, AMP-3D-Wiener achieves more than 2 dB higher PSNR than GPSR, and about 1 dB higher PSNR than TwIST, overall.

Setting 3: In Settings 1 and 2, the measurements are captured with $K=2$ shots. We now test our algorithm on the setting where the number of shots varies from $K=2$ to $K=12$ with pairwise complementary random coded apertures. Specifically, we randomly generate the coded aperture for the $k$-th shot for $k=1,3,5,7,9,11$, and the coded aperture in the $(k+1)$-th shot is the complement of the aperture in the $k$-th shot. In this setting, a moderate amount of noise (20 dB) is added to the measurements. Figure 8 presents the PSNR of the reconstructed cubes as a function of the number of shots, and AMP-3D-Wiener consistently beats GPSR and TwIST.

IV-B Test on natural scenes

Besides the Lego image cube, we have also tested our algorithm on image cubes of natural scenes .The cubes are downloaded from http://personalpages.manchester.ac.uk/staff/ d.h.foster/Hyperspectral$\_$images$\_$of$\_$natural$\_$scenes$\_$04.html and http://per- sonal pages.manchester.ac.uk/staff/d.h.foster/Hyperspectral$\_$images$\_$of$\_$natural$\_$ scenes$\_$02.html. There are two datasets, “natural scenes 2002” and “natural scenes 2004,” each one with 8 image data cubes. The cubes in the first dataset have $L=31$ spectral bands with spatial resolution of around $700\times 700$, whereas the cubes in the second dataset have $L=33$ spectral bands with spatial resolution of around $1000\times 1000$. To satisfy the dyadic constraint of the 2D wavelet, we crop their spatial resolution to be $M=N=512$. Because the spatial dimensions of the cubes “scene 6” and “scene7” in the first dataset are smaller than $512\times 512$, we do not include results for these two cubes.

The measurements are captured with $K=2$ shots, and the measurement rate is $m/n=KM(N+L+1)/(MNL)\approx 0.069$ for “natural scene 2002” and $0.065$ for “natural scene 2004.” We test for measurement noise levels such that the SNRs are 15 dB and 20 dB. The typical runtimes for AMP with 400 iterations, GPSR with 400 iterations, and TwIST with 200 iterations are approximately $2,800$ seconds. The average PSNR over all spectral bands for each reconstructed cube is shown in Tables I and II. We highlight the highest PSNR among AMP-3D-Wiener, GPSR, and TwIST using bold fonts. It can be seen from Tables I and II that AMP-3D-Wiener usually outperforms GPSR by $2-5$ dB in terms of the PSNR, and outperforms TwIST by $0.2-4$ dB, while TwIST outperforms GPSR by up to 3 dB for most of the scenes. Additionally, the results of 6 selected image cubes are displayed in Figure 9 in the form of 2D RGB images.We refer to the tutorial from http://personalpages.manchester.ac.uk/staff/ david.foster/Tutorial$\_$HSI2RGB/Tutorial$\_$HSI2RGB.html and convert 3D image cubes to 2D RGB images. The four rows of images correspond to ground truth, results by AMP-3D-Wiener, results by GPSR, and results by TwIST, respectively. We can see from Figure 9 that the test datasets contain both smooth scenes and scenes with large gradients, and AMP-3D-Wiener consistently reconstructs better than GPSR and TwIST, which suggests that AMP-3D-Wiener is adaptive to various types of scenes.

V Conclusion

In this paper, we considered the compressive hyperspectral imaging reconstruction problem for the coded aperture snapshot spectral imager (CASSI) system. Considering that the CASSI system is a great improvement in terms of imaging quality and acquisition speed over conventional spectral imaging techniques, it is desirable to further improve CASSI by accelerating the 3D image cube reconstruction process. Our proposed AMP-3D-Wiener used an adaptive Wiener filter as a 3D image denoiser within the approximate message passing (AMP) framework. AMP-3D-Wiener was faster than existing image cube reconstruction algorithms, and also achieved better reconstruction quality.

In AMP, the derivative of the image denoiser is required, and the adaptive Wiener filter can be expressed in closed form using a simple formula, and so its derivative is easy to compute. Although the matrix that models the CASSI system is ill-conditioned and may cause AMP to diverge, we helped AMP converge using damping, and reconstructed 3D image cubes successfully. Numerical results showed that AMP-3D-Wiener is robust and fast, and outperforms gradient projection for sparse reconstruction (GPSR) and two-step iterative shrinkage/thresholding (TwIST) even when the regularization parameters for GPSR and TwIST are optimally tuned. Moreover, a significant advantage over GPSR and TwIST is that AMP-3D-Wiener need not tune any parameters, and thus an image cube can be reconstructed by running AMP-3D-Wiener only once, which is critical in real-world scenarios. In contrast, GPSR and TwIST must be run multiple times in order to find the optimal regularization parameters.

Future improvements: In our current AMP-3D-Wiener algorithm for compressive hyperspectral imaging reconstruction, we estimated the noise variance of the noisy image cube within each AMP iteration using (7). In order to denoise the noisy image cube in the sparsifying transform domain, we applied the estimated noise variance value to all wavelet subbands. The noise variance estimation and 3D image denoising method were effective, and helped produce promising reconstruction. However, both the noise variance estimation and the 3D image denoising method may be sub-optimal, because the noisy image cube within each AMP iteration does not contain i.i.d. Gaussian noise, and so the coefficients in the different wavelet subbands may contain different amounts of noise. On the other hand, in the proposed adaptive Wiener filter, the variances of the coefficients in the sparsifying transform domain were estimated empirically within each wavelet subband, whereas it is also possible to apply Wiener filtering via marginal likelihood or generalized cross validation . Therefore, it is possible that the denoising part of the proposed algorithm can be further improved. The study of such denoising methods is left for future work.

In our current AMP-3D-Wiener, the coded apertures must be complementary, because complementary coded apertures ensure that the norm of each column in the matrix ${\bf H}$ in (2) is similar, otherwise, AMP-3D-Wiener may diverge. Although using complementary coded aperture has practical importance, it provides more flexibility in coded aperture design when such a complementary constraint can be removed, and the development for AMP-based algorithms without such constraints is left for future work.

Finally, besides reconstructing image cubes from compressive hyperspectral imaging systems, it would also be interesting to investigate problems such as target detection and unmixing using compressive measurements from hyperspectral imaging systems. We leave these problems for future work.

Acknowledgments

We thank Sundeep Rangan and Phil Schniter for inspiring discussions on approximate message passing; Lawrance Carin, and Xin Yuan for kind help on numerical experiments; Junan Zhu for informative explanations about CASSI systems; Nikhil Krishnan for detailed suggestions on the manuscript; and the reviewers for their careful evaluation of the manuscript.

References