# Compressive Spectral Image Reconstruction using Deep Prior and Low-Rank
Tensor Representation

Jorge Bacca  Yesid Fonseca  and Henry Arguello  Department of Systems
Engineering, Universidad Industrial de Santander, Bucaramanga, Colombia
<EMAIL_ADDRESS>
###### Abstract

Compressive spectral imaging (CSI) has emerged as an alternative spectral
image acquisition technology, which reduces the number of measurements at the
cost of requiring a recovery process. In general, the reconstruction methods
are based on hand-crafted priors used as regularizers in optimization
algorithms or recent deep neural networks employed as an image generator to
learn a non-linear mapping from the low-dimensional compressed measurements to
the image space. However, these data-driven methods need many spectral images
to obtain good performance. In this work, a deep recovery framework for CSI
without training data is presented. The proposed method is based on the fact
that the structure of some deep neural networks and an appropriated low-
dimensional structure are sufficient to impose a structure of the underlying
spectral image from CSI. We analyzed the low-dimension structure via the
Tucker representation, modeled in the first net layer. The proposed scheme is
obtained by minimizing the $\ell_{2}$-norm distance between the compressive
measurements and the predicted measurements, and the desired recovered
spectral image is formed just before the forward operator. Simulated and
experimental results verify the effectiveness of the proposed method.

††journal: osajournal††articletype: Research Article

##  1 Introduction

Spectral imaging (SI) deals with capturing the spatial information of a target
in a broader range of the electromagnetic spectrum compared to a conventional
RGB imaging system. This additional information is useful for some
applications such as biomedical imaging [1], crop identification [2], and
surveillance [3]. SI can be denoted as a 3D tensor
$\bm{\mathcal{X}}\in\mathbb{R}^{M\times N\times L}$ with $M\times N$ as the
spatial pixels and $L$ spectral bands [2]. Traditional methods to acquire SI
are based on scanning along one of its tensor modes, which results in time-
consuming systems, and therefore, prohibits its usage in dynamic scenes [4].

Alternatively, based on the compressive sensing (CS) theory, new imaging
snapshots systems acquire 2D multiplexed projections of a scene instead of
directly acquire all voxels, resulting in an image compression via hardware
[5]. To date, different compressive spectral imaging (CSI) techniques have
been proposed [6, 7, 8, 9, 10, 11]. For instance, the pioneer coded aperture
snapshot spectral imaging (CASSI) system [10] uses optical elements to encode
and disperse the incoming light to acquire 2D intensity projections. Even
though CSI yield efficient sensing, a reconstruction process from the
compressed measurements is needed, since it results in finding a solution to
an under-determined system [5]. This recovery problem is addressed by
representing the 3D scene as a 1D vector and assuming particular spectral
image nature priors in different dimensions used as regularization in an
optimization problem [4, 12]. For instance, [13, 14] assume low total
variation, [9, 7] explore the sparsity assumption of the scene in some
orthogonal basis, [15, 16] use non-local similarity, and [17, 18] employ low-
rank structures. However, these hand-crafted priors do not often represent the
wide variety and non-linearity of spectral images, and the vectorization
ignores the high-dimensional structure of the scene, resulting in low
reconstruction quality [19].

On the other hand, data-driven recovery methods are based on the power of the
deep neural networks as image generators, where the goal is to learn a non-
linear transformation that maps a low-dimensional feature into realistic
spectral images [20]. In particular, with a vast spectral data set, [21, 22,
23, 24] learn inverse networks that map the low-dimensional compressed
measurements to the desired spectral image [25]. These methods have shown high
performance speed and reconstrucion quality. However, they are very dependent
on training data, and small variations in the sensing system would require re-
training of the model [19]. Alternative solutions such as [26], takes the
sensing model into account when solving an optimization problem where the
prior is learned using convolutional auto-encoder with a spectral data set,
and more recently [26, 19, 27, 28] use unrolled-based methods that incorporate
the sensing process into the deep network design, where the prior is
intrinsically learned through end-to-end optimization. Although these methods
have proven to be more general, they still depend on training data.

In this paper, a deep recovery framework for reconstructing spectral images
from CSI measurements without training data requirements is proposed. The
method is based on the fact that the deep convolutional neural networks and
the appropriated low-dimensional representation are sufficient to
learn/generate the image representation without any training data, and
therefore, to recover a spectral image directly from the CSI measurements. In
particular, the proposed method designs a deep neural network, where the first
layer learns a low-dimensional 3D tensor, which is then refined by
convolutional operations to generate the desired reconstruction. The weights
of this neural network are randomly initialized and fitted to guarantee that
the reconstruction suits the CSI measurement via $\ell_{2}$-norm minimization
over the CSI measurement; therefore, the recovered image is formed just before
the forward operator. The proposed method is expressed as an end-to-end
optimization by modeling the forward compressive sensing model as a non-
trainable layer; consequently, it can be solved using any deep learning
algorithm like stochastic gradient descent. Additionally, we analyzed the
importance of the low-dimensional tensor structure in the first layer via low-
rank Tucker representation, which imposes a low-rank 3D-prior. Since there is
no more information available other than the compressive spectral
measurements, the proposed method is more related to hand-crafted techniques.
Results in simulated and real data demonstrate that the proposed method
outperforms the hand-crafted methods in many scenarios and obtains comparable
results with data-driven approaches.

##  2 Related work

###  2.1 Hand-Crafted CS Reconstruction

The traditional CS recovery algorithms are considered hand-designed since they
use some expert knowledge of the signal, known as a signal prior [26]. These
methods are based on optimization techniques that design a data fidelity term,
and incorporate the prior as a regularization term [29]. The most common prior
is assuming that the signal is sparse on a given basis, such as Wavelet [30],
discrete cosine transform (DCT)[5], among others[5]. This sparsity assumption
is imposed in different methods by applying $\ell_{0}$ or $\ell_{1}$
regularizers. Examples of algorithms that use sparsity priors include, the
GPSR [29], ADMM [31], CSALSA[32], ISTA [33], AMP [34] among others. In CSI,
some specific kinds of prior are used. For instance, [9] assumes low total
variation, [7] explores the spatial sparsity assumption of the scene in
Wavelet domain, and the spectral sparsity in the DCT domain [15, 16];
furthermore, [17, 18] employ low-rank structures based on the linear mixture
model. Exploring tensor structure, low-rank tensor recovery methods have been
also proposed [12, 35]. However, these hand-crafted methods require expert
knowledge of the target to select which prior to use. Therefore, they do not
represent the wide variety and the non-linearity of spectral image
representations.

###  2.2 Data-Driven CS Reconstruction

Data-driven recovery methods are based on learning a non-linear inverse
mapping from the compressive measurements to a realistic image. In particular,
with a vast dataset of ground-truth and compressive measurement pairs, these
methods are used to learn a non-linear network by minimizing the distance
between the output of the net and the ground-truth. The main difference
between the state-of-the-art methods is their network architecture. For
instance, [36] learns a stacked auto-encoder, convolution layers are applied
in [37], and convolutional, residual, and fully-connected layers are also used
in [38, 39, 40, 41]. In particular, for CSI, [22] was the first work that used
a data-driven approach, where, an initialization obtained from TwiST [42] was
refined using denoising networks; [19] proposed a particular model to explore
the spatial and spectral information and to design the coded aperture usually
included in CSI architectures. Furthermore, based on the structure of the
U-net, [24] proposed a non-linear mapping replacing the 2D for 3D
convolutions, and [23] developed a generative model based on the U-net. These
methods have shown high performance in reconstruction quality, and once
trained, they allow real-time reconstruction. However, these approaches are
highly dependent on the data-set used. Furthermore, small-variations in the
compressive measurements, such as type of noise or changes in the sensing
matrix, would require a time-consuming re-training.

###  2.3 Non-Linear Image Priors for CS Methods

Recently, some works have considered the sensing model to proposed a mixed
approach which considers the hand crafted as well as the data-driven CS
reconstruction. In particular, these methods use a deep network or denoiser to
replace the hand-crafted prior, then, this non-linear prior is employed in the
optimization algorithm [38]. For instance, Plug-and-play priors (PnP) use pre-
existing denoisers as a proximal step [43, 44], [45] learns the proximal
mapping using a convolutional network, and [26] learns a SI prior, through a
convolutional autoencoder, which is then incorporated into the optimization
problem. More recently, D-AMP [46], ISTA-Net [47], ADMM-Net [48], and DNU [28]
use the unrolled based method that incorporates the optimization steps into
the deep network architecture using residual networks; consequently, they can
learn the prior and the parameters via end-to-end training. This strategy is
also employed for CSI in [19, 27]. Although these methods have proven to be
more general, they still depend on training data, which is limited in SI.

###  2.4 Deep Image Prior using Generative Model

The generative model (GM) has been used for CS recovery [49]. The goal in GM
is to generate a realistic image from a low-dimensional latent representation.
For instance, [49, 50] use a pre-trained deep neural network and obtains the
low-dimensional representation, which minimizes the distance between the
compressive measurements and the output of the net. On the other hand, [51]
shows that a pre-trained network is not necessary. Instead of finding the low-
dimensional latent space, [51] uses a fixed random variable as latent space,
then the weights of the model are updated to obtain an optimal result. The
drawback of this method is its sensitivity to changes in the application, the
fixed latent space or the network architecture, which usually require small
random disturbances to obtain a good performance. The proposed method in this
work is closely related to [50, 51], where the parameters of the network model
and the low-dimensional representation (based on a Tucker representation,
which is useful for SI) are optimized in an end-to-end approach for a CSI
architecture.

### Notation:

Through the paper, vectors are represented with boldface lowercase letters,
e.g., $\bm{x},$ and matrices are denoted as boldface capital letters
$\mathbf{X}$. The 3D tensors are denoted as
$\bm{\mathcal{X}}\in\mathbb{R}^{M\times N\times L}$ and the $1$-mode product
of a tensor $\bm{\mathcal{X}}_{o}\in\mathbb{R}^{M_{p}\times N_{p}\times
L_{p}}$ with a matrix $\mathbf{U}\in\mathbb{R}^{M\times M_{p}}$ is written as
$\bm{\mathcal{X}}=\bm{\mathcal{X}}_{o}\times_{1}\mathbf{U}$ where
$\bm{\mathcal{X}}\in\mathbb{R}^{M\times N_{p}\times L_{p}}$, and

$\displaystyle\bm{\mathcal{X}}_{(m,n,\ell)}=\sum_{\hat{m}=1}^{M_{p}}\mathbf{U}_{(m,\hat{m})}\bm{\mathcal{X}}_{o(\hat{m},n,\ell)}.$

In the same way, the 2-mode and 3-mode products can be defined. We introduce
the function $\text{shift}_{\ell}(\cdot):\mathbb{R}^{M\times
N}\rightarrow\mathbb{R}^{M\times(N+L-1)}$ which refers to a shifting operator,
i.e., for a given $\mathbf{X}$ we have that

$\displaystyle\text{shift}_{\ell}(\mathbf{X}):=\begin{cases}\mathbf{X}_{(m,n-\ell)},&\text{
if }1\leq n-l\leq N\\\ 0,&\text{ otherwise}.\end{cases}$

Finally, the function $\text{vect}(\cdot):\mathbb{R}^{M\times N\times
L}\rightarrow\mathbb{R}^{MNL}$ represents the vectorization of a tensor.

##  3 Compressed Measurements Acquisition

Figure 1: Physical sensing phenomena in CASSI, which is the CSI prototype used
to validate the proposed approach.

The CASSI sensing approach is used in order to acquire the compressed
measurements of a spectral scene [10]. This architecture is composed of three
main optical elements: a coded aperture, a prism as a dispersive element, and
a gray-scale detector, as illustrated in Fig 1. The spatial-spectral data cube
is represented as $\bm{\mathcal{X}}\in\mathbb{R}^{M\times N\times L}$ with
$M\times N$ spatial dimensions, $L$ spectral bands, and
$\mathbf{X}_{\ell}\in\mathbb{R}^{M\times N}$ denotes the 2D spectral intensity
image of $\bm{\mathcal{X}}$ at the $\ell$-th spectral band. As shown in Fig.
1, each spatial position of the scene is modulated by a coded aperture
$\mathbf{C}\in\\{0,1\\}^{M\times N}$, which block/unblock the incoming light,
then, the coded spectral scene passes through the prism creating a horizontal
shifting. Finally, the coded shifted spectral scene is integrated along the
spectral axis by the detector, resulting in the 2D compressed measurement
$\mathbf{Y}\in\mathbb{R}^{M\times(N+L-1)}$. In CSI, it is possible to acquire
$S<L$ different measurement snapshots of the same spectral data cube employing
$S$ different patterns in the coded aperture. Therefore, the output of the
sensing process at the $s$-th spectral snapshot can be mathematically
expressed as

$\mathbf{Y}^{(s)}=\sum_{\ell=1}^{L}\text{shift}_{\ell-1}\left(\mathbf{X}_{\ell}\odot\mathbf{C}^{(s)}\right),$
(1)

where the $\ell$-th spectral band, $\mathbf{X}_{\ell}$, of the tensor
$\bm{\mathcal{X}}$ is shifted with the operator
$\text{shift}_{\ell-1}(\cdot)$, and $\odot$ denotes the element-wise product
with the 2D coded aperture $\mathbf{C}^{(s)}$.

The CASSI sensing model can be seen as a linear operator, after stacking the
measurements of multiple shots as
$\bm{y}=[\text{vect}(\mathbf{Y}^{(1)})^{T},\cdots\text{vect}(\mathbf{Y}^{(S)})^{T}]$.
Thus, the system matrix model can be expressed as

$\bm{y}=\mathbf{H}\text{vect}(\bm{\mathcal{X}}),$ (2)

where $\mathbf{H}\in\mathbb{R}^{SM(N+L-1)\times MNL}$ represents the linear
sensing matrix of CASSI.

##  4 Compressive Spectral Reconstruction

Figure 2: Visual representation of the proposed deep neural scheme, where the
boxes with background color represent the learning parameters, the white box
stand for the non-trainable CSI system, and the non-box blocks represent the
outputs of the layers.

The goal in CSI is to recover the spectral image
$\bm{\mathcal{X}}\in\mathbb{R}^{M\times N\times L}$ from the compressive
measurements $\bm{y}$. Since $SM(N+L-1)\ll MNL$, this problem consists in
solving an undetermined system, which is addressed by restricting the feasible
set of solutions using image priors as regularizers. A tensor formulation for
addressing this problem is described below

$\displaystyle\underset{\bm{\mathcal{Z}}_{o}\in\mathbb{R}^{M\times N\times
L}}{\mbox{minimize}}\hskip 8.53581pt$
$\displaystyle\frac{1}{2}\left\|\bm{y}-\mathbf{H}\text{vect}\left(\bm{\mathcal{X}}\right)\right\|^{2}_{2}+\lambda\cdot\phi(\bm{\mathcal{Z}^{\prime}}_{o})$
(3) subject to
$\displaystyle\bm{\mathcal{X}}=\bm{\mathcal{Z}^{\prime}}_{o}\times_{1}\mathbf{U}^{\prime}\times_{2}\mathbf{V}^{\prime}\times_{3}\mathbf{W}^{\prime},$

where the matrices $\mathbf{U}^{\prime}\in\mathbb{R}^{M\times
M},\mathbf{V}^{\prime}\in\mathbb{R}^{N\times N}$ and
$\mathbf{W}^{\prime}\in\mathbb{R}^{L\times L}$ are fixed and known orthogonal
matrices, which usually are the matrix representation of the Wavelet and the
Discrete Cosine transforms; $\bm{\mathcal{Z}^{\prime}}_{o}$ is the
representation of the spectral image in the given basis and
$\phi(\cdot):\mathbb{R}^{M\times N\times L}\rightarrow\mathbb{R}$ is a
regularization function that imposes particular image priors with $\lambda$ as
the regularization parameter [29].

Unlike the hand-craft priors as sparsity [5], we explore the power of some
deep neural networks as image generators that map a low-dimensional feature
tensor $\bm{\mathcal{Z}}\in\mathbb{R}^{M\times N\times L}$ to the image as

$\bm{\mathcal{X}}=\mathcal{M}_{\bm{\theta}}(\bm{\mathcal{Z}}),$ (4)

where $\mathcal{M}_{\bm{\theta}}(\cdot)$ represents a deep network, with
${\bm{\theta}}$ as the net-parameters. To ensure a low-dimensional structure
over the feature tensor, this work used the Tucker representation, i.e.,
$\bm{\mathcal{Z}}=\bm{\mathcal{Z}}_{o}\times_{1}\mathbf{U}\times_{2}\mathbf{V}\times_{3}\mathbf{W}$
with $\bm{\mathcal{Z}}_{o}\in\mathbb{R}^{M_{\rho}\times N_{\rho}\times
L_{\rho}}$ as a 3D low dimensional tensor, with $M_{\rho}<M$, $N_{\rho}<N$ and
$L_{\rho}<L$. This representation maintains the 3D structure of the spectral
images, exploits the inherent low-rank of this data [52, 53], and also
implicitly constraint the output $\bm{\mathcal{X}}$ in a low-dimensional
manifold via the architecture and the weights of the net [50].

In this paper, we are focused in a blind representation, where instead of have
a pre-training network or huge amount of data to train this deep neural
representation, we express an optimization problem which learns the weight
$\bm{\theta}$ in the generative network $\mathcal{M}_{\bm{\theta}}$ and also
the tensor feature $\bm{\mathcal{Z}}$ with its Tucker representation elements
as $\bm{\mathcal{Z}}_{o},\mathbf{U},\mathbf{V}$ and $\mathbf{W}$. All the
parameters of this optimization problem are randomly initialized and the only
available information are the compressive measurements and the sensing model,
i.e, the optimization problem is data training independent. In particular, we
explore the prior implicitly captured by the choice of the generator network
structure, which is usually composed of convolutional operations, and the
importance of the low-rank representation feature, therefore, the proposed
method consists of solving the following optimization problem

$\displaystyle\underset{{\bm{\theta}},\bm{\mathcal{Z}}_{o},\mathbf{U,V,W}}{\mbox{minimize}}\hskip
8.53581pt$
$\displaystyle\frac{1}{2}\left\|\bm{y}-\mathbf{H}\text{vect}\left(\mathcal{M}_{{\bm{\theta}}}\left(\bm{\mathcal{Z}}\right)\right)\right\|^{2}_{2}$
(5) subject to
$\displaystyle\bm{\mathcal{Z}}=\bm{\mathcal{Z}}_{o}\times_{1}\mathbf{U}\times_{2}\mathbf{V}\times_{3}\mathbf{W},$

where the recovery is
$\bm{\mathcal{X}}^{*}=\mathcal{M}_{\bm{\theta}^{*}}(\bm{\mathcal{Z}}_{o}^{*}\times_{1}\mathbf{U}^{*}\times_{2}\mathbf{V}^{*}\times_{3}\mathbf{W}^{*})$.
This optimization problem can be solved using an end-to-end neural network
framework, as shown in Fig. 2. In this way, the input, that is common in all
neural networks, is replaced with a custom layer with
$\bm{\mathcal{Z}}_{o},\mathbf{U,V,W}$ as learnable parameters, which construct
the low-rank Tucker representation of $\bm{\mathcal{Z}}$, then this tensor
$\bm{\mathcal{Z}}$ is refined with convolutional layers via
$\mathcal{M}_{\bm{\theta}}(\bm{\mathcal{Z}})$; these optimization variables
are represented by the first two blue-blocks in the Fig. 2. The final layer in
the proposed method is a non-training layer which models the forward sensing
operator
$\mathbf{H}\text{vect}\left(\mathcal{M}_{{\bm{\theta}}}\left(\bm{\mathcal{Z}}\right)\right)$
to obtain the compressive measurements $\bm{y}$ as the output of the net.
Therefore, the problem in (5) can be solved with state-of-the-art deep
learning optimization algorithm, such as, stochastic gradient descent. Once
the parameters are optimized, the desired SI is recovered just before the non-
trainable layer labeled as "CSI system" in Fig. 2.

##  5 Simulation and Results

In this section, the performance of the proposed compressive spectral image
reconstruction approach is presented. The performance metrics used are the
peak-signal-to-noise (PSNR) [5], the structural similarity (SSIM) [54], and
the spectral angle mapping (SAM) [17]. PSNR and SSIM are calculated as the
average of each 2D spatial image through the bands, and the SAM is the average
of all spectral pixels. Four different tests are presented to validate the
proposed method. The first test evaluates the importance of the low-rank
tensor representation; the second test compares the recovery of the data-
driven methods with the proposed method; the third evaluates the proposed
method in different noisy scenarios and for a different number of shots
against the non-data dependent state-of-the-art algorithms; and finally, the
proposed method is evaluated using two compressive spectral images obtained
with a real test-bed implementation 111The code can be find
https://github.com/jorgebaccauis/Deep_Prior_Low_Rank.

###  5.1 Rank level

Figure 3: Visual representation of the three network models used: U-Net-based,
AutoencoderNet-based and ResNet-based. The color represents the different
layers in each network.  Figure 4: PSNR Box plot for the different network
architectures varying the rank factor $\rho$, with $5$ run trials.

This section evaluates the importance of the rank level in the 3D tensor using
the Tucker representation, which is placed at the first block of our model, as
illustrated in Fig. 2. For that, two spectral images with $M=256\times N=256$
pixels, and $L=10$ spectral bands between $400$ and $700$nm from [55] where
chosen. Three different network architectures were tested as ”Convolutional
Layers” for the second block in the Figure 2. The first network architecture
is a simple ResNet-based model [56], with a single skip connection and four
convolutional layers, as shown in the Figure 3 with $2.150$ parameters. The
second architecture, also shown in Fig. 3, is a convolution Autoencoder-based
[57], with $8.160$ training parameters, and six convolutional layers. The
third architecture tested and depicted in FIg. 3, is a Unet-based [58],
without drop-out layers, and, in the contracting part, the feature information
is increased using multiples of $L=10$, i.e., $L,2L$ and $3L$ as is
illustrated in Fig.3, resulting in $92.190$ training parameters. This test is
focused on a single snapshot for a randomly coded aperture generated from a
Bernoulli distribution with mean 0.5.

As mentioned, the tensor feature $\bm{\mathcal{Z}}\in\mathbb{R}^{M\times
N\times L}$ comes from a low-dimensional kernel
$\bm{\mathcal{Z}}_{o}\in\mathbb{R}^{M_{p}\times N_{p}\times L_{p}}$; then, to
evaluate the importance of the rank-level in the Tucker representation, we
establish the following relationship

$\frac{M_{p}}{M}=\frac{N_{p}}{N}=\frac{L_{p}}{L}=\rho,$ (6)

where $\rho\in(0,1]$, is referred as the hyper-parameter rank factor.
Furthermore, as the parameters of the problem in (5) are randomly initialized,
we simulated five realizations. The average results for this 5 realizations
are summarized in the Figure 4. Notice that for the three network
architectures and the two datasets, the rank factor is a crucial hyper-
parameter to obtain a good reconstruction. In particular, for the DataSet 1
the optimal value is $\rho=\\{0.5,0.7,0.6\\}$ for the ResNet-based,
AutoeconderNet-based, and Unet-based, respectively; and
$\rho=\\{0.4,0.6,0.7\\}$ for DataSet 2. Furthermore, notice that a small value
of $\rho$ presents the worst case for all the networks. Also, notice all the
network configurations obtain around 31 dB, which is the best-obtained
results, for different $\rho$ values; however, the AutoencoderNet-based is
more stable compared with the other networks. This result shows the importance
of the low-rank tensor representation in the first layer, where the optimal
value changes for each dataset and each network architecture.

###  5.2 Data-Driven Methods Comparison

Figure 5:  Two reconstructed scenes using the 5 learning-based methods and the
two variations of the proposed method, i.e., (AutoEncoder, UNet)-Based.

Although the proposed method does not need data to work, this test compares
its results with the data-driven approaches to demonstrate the quality
achieved. In particular, we use five learning-based methods for comparison:
HSCNN [22], ISTA-Net [47], Autoencoder [26]; HIR-DSSP [19] and DNU [28]. These
methods where trained using the public ICVL [59], Harvard [60], and KAIST [26]
hyperspectral image data-sets using their available codes; the sensing process
was evaluated for a single snapshot, according to [28]. For the proposed
method, the two network architectures were evaluated, i.e., AutoEnconder-
Based, and UNet-based. Two testing images of $512\times 512$ of spatial
resolution and $31$ spectral bands were chosen to evaluate the different
methods, and the reconstruction results and ground truth are shown in Fig. 5.
It can be observed that the two variants of the proposed method outperform in
visual and quantitative results to HSCNN, ISTA-Net, AutoEnconder, HIR-DSSP, up
to $(5/0.030/0.020)$ in terms of (PSNR/SSIM/SAM),respectively, and show
comparable/close results with respect to the DNU method, which is the best
data-driven method. However, the proposed method has the advantage that it
does not require training data compared with the driven-data methods, i.e.,
only the compressive measurements are available for the proposed approach.

Table 1: Mean performance comparison for the different recovery methods
varying the number of snapshots and noise in SNR dB.

Shots | Noise | Metrics | GPSR | ADMM | CSALSA |  | PnP
---
ADMM
DIP | Prop.
1 | $\infty$ | PSNR | 25.66 | 24.32 | 25.59 | 28.99 | 27.93 | 30.92
SSIM | 0.701 | 0.726 | 0.790 | 0.860 | 0.766 | 0.874
SAM | 0.145 | 0.108 | 0.152 | 0.060 | 0.089 | 0.055
30 | PSNR | 25.52 | 22.68 | 25.46 | 28.82 | 27.19 | 29.29
SSIM | 0.699 | 0.653 | 0.701 | 0.844 | 0.772 | 0.864
SAM | 0.156 | 0.112 | 0.167 | 0.073 | 0.089 | 0.062
20 | PSNR | 24.67 | 21.45 | 22.19 | 25.42 | 27.53 | 27.94
SSIM | 0.682 | 0.625 | 0.672 | 0.713 | 0.783 | 0.794
SAM | 0.210 | 0.220 | 0.195 | 0.138 | 0.084 | 0.080
2 | $\infty$ | PSNR | 28.15 | 27.14 | 29.45 | 32.28 | 30.578 | 32.88
SSIM | 0.802 | 0.786 | 0.823 | 0.916 | 0.867 | 0.918
SAM | 0.095 | 0.120 | 0.098 | 0.041 | 0.068 | 0.039
30 | PSNR | 27.95 | 26.15 | 29.02 | 31.68 | 30.39 | 31.99
SSIM | 0.789 | 0.726 | 0.835 | 0.903 | 0.878 | 0.910
SAM | 0.099 | 0.130 | 0.099 | 0.055 | 0.066 | 0.053
20 | PSNR | 26.34 | 23.41 | 24.49 | 27.48 | 30.06 | 30.56
SSIM | 0.752 | 0.687 | 0.701 | 0.764 | 0.865 | 0.880
SAM | 0.143 | 0.247 | 0.201 | 0.125 | 0.071 | 0.070
3 | $\infty$ | PSNR | 29.99 | 28.25 | 30.80 | 34.40 | 32.21 | 34.30
SSIM | 0.820 | 0.775 | 0.880 | 0.949 | 0.894 | 0.950
SAM | 0.092 | 0.098 | 0.075 | 0.031 | 0.064 | 0.031
30 | PSNR | 29.53 | 27.50 | 30.78 | 33.19 | 31.54 | 33.45
SSIM | 0.801 | 0.758 | 0.875 | 0.925 | 0.910 | 0.930
SAM | 0.090 | 0.101 | 0.089 | 0.047 | 0.057 | 0.040
20 | PSNR | 27.98 | 24.42 | 25.26 | 28.10 | 29.27 | 31.15
SSIM | 0.775 | 0.685 | 0.695 | 0.778 | 0.883 | 0.900
SAM | 0.114 | 0.220 | 0.205 | 0.124 | 0.066 | 0.062
4 | $\infty$ | PSNR | 30.59 | 30.78 | 31.24 | 38.15 | 32.60 | 36.42
SSIM | 0.889 | 0.891 | 0.921 | 0.973 | 0.914 | 0.966
SAM | 0.072 | 0.062 | 0.043 | 0.022 | 0.057 | 0.29
30 | PSNR | 30.02 | 28.88 | 30.79 | 34.44 | 32.57 | 34.21
SSIM | 0.832 | 0.804 | 0.820 | 0.939 | 0.847 | 0.939
SAM | 0.085 | 0.902 | 0.062 | 0.043 | 0.076 | 0.042
20 | PSNR | 28.16 | 25.67 | 26.74 | 28.33 | 32.30 | 33.74
SSIM | 0.798 | 0.721 | 0.742 | 0.776 | 0.909 | 0.925
SAM | 0.130 | 0.175 | 0.201 | 0.127 | 0.059 | 0.102

###  5.3 Robustness Analysis

Figure 6:  Two RGB false color reconstructed scenes using the non-data driven
methods and the proposed method with its respective metrics are presented.
Additionally, the ground-truth and a spectral point of each scene is shown.

Numerical simulations were conducted to demonstrate the robustness of the
proposed method at different levels of additive Gaussian noise and the number
of shots, using the two spectral image obtained in [55]. It is well known that
the data-driven methods distribution of training and test data must be similar
to obtain good results, for this reason, in this experiment, the proposed
method was compared with the state-of-art non-data driven methods.
Specifically, we compare the proposed method with the GPSR [29], using
sparsity assumption in the Wavelet Kronecker Discrete Cosine transform
implemented as in [8], ADMM [31] using low-rank prior implemented as in [17],
CSALSA [32] using the 3D total variation, PnP-ADMM [43] using BM3D as
denoiser, and Deep Image Prior [51] using the ResNet based network. Three
different noise levels were evaluated: 20, 30 dB of signal to noise ratio
(SNR) and noiseless case that results in $\infty$ dB. Further, the number of
snapshots is varied between 1, 2, 3, and 4 snapshots using the CASSI system,
expressed mathematically as in (2). For this experiment, the ResNet-based
network was used as the "Convolutional layers" in the proposed model, and the
rank factor was fixed as $\rho=0.5$ and $\rho=0.4$ for the DataSet 1 and
DataSet 2, respectively. Table LABEL:table:Robust_analysis, presents a
comparison of the performance in terms of PSNR, SSIM, and SAM metrics, for the
different methods (The results are the average of the two DataSet). Boldface
indicates the best result for each case, and the second-best result is
underlined. From the Table LABEL:table:Robust_analysis, it can be seen that
the proposed method outperforms in almost all cases the other methods.
Furthermore, the proposed method shows better noise robustness compared to the
other methods, since the maximum quality difference between the noise levels
studied is 3 dB, compared to 3dB, 5dB, 5dB, 5dB, 10dB and 3dB for the GPSR,
ADMM, CSALSA, PnP -ADMM and Deep Image Prior (DIP), respectively.
Additionally, as expected, when the number of snapshots per image increases,
all the methods improve their reconstruction quality; in particular, the
difference between 1 and 4 snapshots for the noiseless case, is up to, 5dB,
6dB, 6dB, 10dB 5dB and 6dB for GPSR, ADMM, CSALSA, PnP-ADMM, DIP and the
proposed method, respectively.

To visualize the reconstructions and analyze the results in more detail,
Figure 6 shows an RGB false color for the reconstruction of each method, for a
single CASSI shot, which is the extreme case in terms of compression. Note,
that the proposed method, in the zoomed insets, is much cleaner than its
counterparts. Additionally, to see the behavior, a single spatial point of
each reconstruction for the two Datasets are also presented in Figure 6. It
can be seen that the spectral signatures obtained by the proposed method
closely resemble the ground-truth.

###  5.4 Validation in a Real Testbed Implementation

This section evaluates the proposed method with real measurements acquired
using a testbed implementation. For this section, the ResNet-based model was
used with ($\rho=0.4$), and learning rate $1e-3$. Specifically, two different
scenarios of compressed projections were assessed, which are described as
follows.

Figure 7: Testbed CASSI implementation where the relay lens focuses the
encoded light by the DMD into the sensor after dispersed by the prism.

####  5.4.1 Binary Coded Aperture

This scenario was carried out for one snapshot of the CASSI testbed laboratory
implementation depicted in Fig. 7. This setup contains a $100$-$nm$ objective
lens, a high-speed digital micro-mirror device (DMD) (Texas Instruments-
DLI4130), with a pixel size of $13,6\mu m$, where the CA is implemented, an
Amici Prism (Shanghai Optics), and a CCD (AVT Stingray F-145B) camera with
spatial resolution $1388\times 1038$, and pitch size of $6.45\mu m$. The CA
spatial distribution for the snapshot comes from blue noise patterns, i.e.,
this CA is designed according to [61]. Notice that the robustness analysis
summarized in Table LABEL:table:Robust_analysis, showed that the three best
recovery methods were the PnP-ADMM, DIP, and the proposed method; therefore,
we decided also to compare them using this real data.

Figure 8:  (Left) RGB visual representation of the scene obtained with the
different methods, (Right), two spectral signatures of the recovered scenes.
Figure 9:  (Top) RGB visual representation of the scene obtained with the GPSR
method used in [62] and the proposed method, (Bottom), normalized spectral
signatures of the recovered scenes.

Figure 8 presents the RGB scene obtained with a traditional camera, and the
false-colored RGB images corresponding to reconstructed spectral images using
the different solvers. Furthermore, the spectral responses of two particular
spatial locations in the scene, indicated as red points in the images, are
also included and compared with the spectral behavior using a commercially
available spectrometer (Ocean Optics USB2000+). The visual results show that
the proposed method yield better spatial and spectral reconstruction since the
RGB reconstructed is sharper in the proposed scheme, and the spectral
signatures are closer to those taken by the spectrometer, this is, the SAM of
the normalized signatures obtained from the PnP-ADMM algorithm is 0.188, Deep
Image Prior is 0.205, and the SAM associated to the proposed method is 0.120.
These numerical results validate the performance of the proposed method with
real data for a real CASSI setup using a binary-coded aperture.

####  5.4.2 Colored Coded Aperture

The real data for this second test was provided by [62]. In particular, the
main difference with the data of Section 5.4.1 is that the spatial modulation
is a Colored CA, where each pixel can be seen as a filter with its spectral
response, (further details regarding Colored CA can be found in [8, 62]). The
optical elements in this testbed implementation were the same used in the
previous setup, where the DMD was used to emulate the Colored CA. The coding
and the scene were implemented to have a spatial resolution of $256\times 256$
pixels and $L=8$ as the resolvable bands. [62] uses a hand-crafted method,
which does not require training data and the GPSR algorithm was used as a
recovery algorithm; therefore, the proposed method was compared with this
method. Figure 9(Top) shows the RGB mapping of the recovered scene. There, it
can be seen that the proposed method provides a cleaner version of the scene.
Additionally, two spatial points were chosen to evaluate the spectral
behavior. The normalized spectral response of the proposed method and the GPSR
algorithm are also included in the Figure 9(Bottom). It can be seen that the
spectral signature provided by the proposed method is closer to the obtained
with the spectrometer compared with the GPSR method, in fact, the SAM of the
normalized signatures obtained from the GPSR algorithm is 0.120 and the SAM
associated to the proposed method is 0.057. These results validate the
effectiveness of the proposed method on real data for two variations of CASSI
systems.

##  6 Conclusions

A method for reconstructing spectral images from the CSI measurements has been
proposed. The proposed scheme is based on the fact that the spectral images
can be generated from a convolutional network whose input features comes from
a low-rank Tucker representation. Although the proposed method is based on a
convolutional network framework, it does not require training data, only the
compressed measurements. This method was evaluated in three scenarios:
noiseless, noisy, and real data implementation. In all of them, the proposed
method outperforms the image quality reconstruction compared with state-of-
the-art methods. In particular, the proposed method with 20 SNR levels of
noise in the CSI measurements outperforms its counterparts in up to 4 dB in
the PSNR measure. Although the proposed method was tested in two variation of
CSI system, it can be extended and used in compressive systems where the data
set is limited.

## References

  * [1] G. Lu and B. Fei, “Medical hyperspectral imaging: a review,” Journal of biomedical optics 19, 010901 (2014).
  * [2] L. Zhang, L. Zhang, D. Tao, and X. Huang, “Tensor discriminative locality alignment for hyperspectral image spectral–spatial feature extraction,” IEEE Transactions on Geoscience and Remote Sensing 51, 242–256 (2012).
  * [3] P. W. Yuen and M. Richardson, “An introduction to hyperspectral imaging and its application for security, surveillance and target acquisition,” The Imaging Science Journal 58, 241–253 (2010).
  * [4] C. Hinojosa, J. Bacca, and H. Arguello, “Coded aperture design for compressive spectral subspace clustering,” IEEE Journal of Selected Topics in Signal Processing 12, 1589–1600 (2018).
  * [5] G. R. Arce, D. J. Brady, L. Carin, H. Arguello, and D. S. Kittle, “Compressive coded aperture spectral imaging: An introduction,” IEEE Signal Processing Magazine 31, 105–115 (2014).
  * [6] X. Cao, T. Yue, X. Lin, S. Lin, X. Yuan, Q. Dai, L. Carin, and D. J. Brady, “Computational snapshot multispectral cameras: Toward dynamic capture of the spectral world,” IEEE Signal Processing Magazine 33, 95–108 (2016).
  * [7] C. V. Correa, C. Hinojosa, G. R. Arce, and H. Arguello, “Multiple snapshot colored compressive spectral imager,” Optical Engineering 56, 041309 (2016).
  * [8] H. Arguello and G. R. Arce, “Colored coded aperture design by concentration of measure in compressive spectral imaging,” IEEE Transactions on Image Processing 23, 1896–1908 (2014).
  * [9] A. Wagadarikar, R. John, R. Willett, and D. Brady, “Single disperser design for coded aperture snapshot spectral imaging,” Applied optics 47, B44–B51 (2008).
  * [10] M. Gehm, R. John, D. Brady, R. Willett, and T. Schulz, “Single-shot compressive spectral imaging with a dual-disperser architecture,” Optics express 15, 14013–14027 (2007).
  * [11] S. Shauli, O. Yaniv, A. Marwan, A. Ibrahim, D. G. Blumberg, and A. Stern, “Dual-camera design for hyperspectral and panchromatic imaging, using a wedge shaped liquid crystal as a spectral multiplexer,” Scientific Reports (Nature Publisher Group) 10 (2020).
  * [12] S. Zhang, L. Wang, Y. Fu, X. Zhong, and H. Huang, “Computational hyperspectral imaging based on dimension-discriminative low-rank tensor recovery,” in _Proceedings of the IEEE International Conference on Computer Vision,_ (2019), pp. 10183–10192.
  * [13] D. Kittle, K. Choi, A. Wagadarikar, and D. J. Brady, “Multiframe image estimation for coded aperture snapshot spectral imagers,” Applied Optics 49, 6824–6833 (2010).
  * [14] L. Wang, Z. Xiong, D. Gao, G. Shi, and F. Wu, “Dual-camera design for coded aperture snapshot spectral imaging,” Applied optics 54, 848–858 (2015).
  * [15] Y. Fu, Y. Zheng, I. Sato, and Y. Sato, “Exploiting spectral-spatial correlation for coded hyperspectral image restoration,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition,_ (2016), pp. 3727–3736.
  * [16] L. Wang, Z. Xiong, G. Shi, F. Wu, and W. Zeng, “Adaptive nonlocal sparse representation for dual-camera compressive hyperspectral imaging,” IEEE transactions on pattern analysis and machine intelligence 39, 2104–2111 (2016).
  * [17] J. Bacca, C. V. Correa, and H. Arguello, “Noniterative hyperspectral image reconstruction from compressive fused measurements,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 12, 1231–1239 (2019).
  * [18] T. Gelvez, H. Rueda, and H. Arguello, “Joint sparse and low rank recovery algorithm for compressive hyperspectral imaging,” Applied optics 56, 6785–6795 (2017).
  * [19] L. Wang, C. Sun, Y. Fu, M. H. Kim, and H. Huang, “Hyperspectral image reconstruction using a deep spatial-spectral prior,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition,_ (2019), pp. 8032–8041.
  * [20] R. Hyder and M. S. Asif, “Generative models for low-rank video representation and reconstruction from compressive measurements,” in _2019 IEEE 29th International Workshop on Machine Learning for Signal Processing (MLSP),_ (IEEE, 2019), pp. 1–6.
  * [21] L. Wang, T. Zhang, Y. Fu, and H. Huang, “Hyperreconnet: Joint coded aperture optimization and image reconstruction for compressive hyperspectral imaging,” IEEE Transactions on Image Processing 28, 2257–2270 (2018).
  * [22] Z. Xiong, Z. Shi, H. Li, L. Wang, D. Liu, and F. Wu, “Hscnn: Cnn-based hyperspectral image recovery from spectrally undersampled projections,” in _Proceedings of the IEEE International Conference on Computer Vision Workshops,_ (2017), pp. 518–525.
  * [23] X. Miao, X. Yuan, Y. Pu, and V. Athitsos, “$\lambda$-net: Reconstruct hyperspectral images from a snapshot measurement,” in _IEEE/CVF Conference on Computer Vision (ICCV),_ vol. 1 (2019).
  * [24] D. Gedalin, Y. Oiknine, and A. Stern, “Deepcubenet: reconstruction of spectrally compressive sensed hyperspectral images with deep neural networks,” Optics Express 27, 35811–35822 (2019).
  * [25] J. Bacca, L. Galvis, and H. Arguello, “Coupled deep learning coded aperture design for compressive image classification,” Optics Express 28, 8528–8540 (2020).
  * [26] I. Choi, D. S. Jeon, G. Nam, D. Gutierrez, and M. H. Kim, “High-quality hyperspectral reconstruction using a spectral prior,” ACM Transactions on Graphics (TOG) 36, 1–13 (2017).
  * [27] T. Zhang, Y. Fu, L. Wang, and H. Huang, “Hyperspectral image reconstruction using deep external and internal learning,” in _Proceedings of the IEEE International Conference on Computer Vision,_ (2019), pp. 8559–8568.
  * [28] L. Wang, C. Sun, M. Zhang, Y. Fu, and H. Huang, “Dnu: Deep non-local unrolling for computational spectral imaging,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition,_ (2020), pp. 1661–1671.
  * [29] M. A. Figueiredo, R. D. Nowak, and S. J. Wright, “Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems,” IEEE Journal of selected topics in signal processing 1, 586–597 (2007).
  * [30] E. J. Candès and M. B. Wakin, “An introduction to compressive sampling,” IEEE signal processing magazine 25, 21–30 (2008).
  * [31] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein _et al._ , “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends in Machine learning 3, 1–122 (2011).
  * [32] M. V. Afonso, J. M. Bioucas-Dias, and M. A. Figueiredo, “An augmented lagrangian approach to the constrained optimization formulation of imaging inverse problems,” IEEE Transactions on Image Processing 20, 681–695 (2010).
  * [33] I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences 57, 1413–1457 (2004).
  * [34] D. L. Donoho, A. Maleki, and A. Montanari, “Message-passing algorithms for compressed sensing,” Proceedings of the National Academy of Sciences 106, 18914–18919 (2009).
  * [35] S. Yang, M. Wang, P. Li, L. Jin, B. Wu, and L. Jiao, “Compressive hyperspectral imaging via sparse tensor and nonlinear compressed sensing,” IEEE Transactions on Geoscience and Remote Sensing 53, 5943–5957 (2015).
  * [36] A. Mousavi, A. B. Patel, and R. G. Baraniuk, “A deep learning approach to structured signal recovery,” in _2015 53rd annual allerton conference on communication, control, and computing (Allerton),_ (IEEE, 2015), pp. 1336–1343.
  * [37] A. Mousavi and R. G. Baraniuk, “Learning to invert: Signal recovery via deep convolutional networks,” in _2017 IEEE international conference on acoustics, speech and signal processing (ICASSP),_ (IEEE, 2017), pp. 2272–2276.
  * [38] A. Dave, A. K. Vadathya, R. Subramanyam, R. Baburajan, and K. Mitra, “Solving inverse computational imaging problems using deep pixel-level prior,” IEEE Transactions on Computational Imaging 5, 37–51 (2018).
  * [39] H. Palangi, R. Ward, and L. Deng, “Distributed compressive sensing: A deep learning approach,” IEEE Transactions on Signal Processing 64, 4504–4518 (2016).
  * [40] H. Yao, F. Dai, S. Zhang, Y. Zhang, Q. Tian, and C. Xu, “Dr2-net: Deep residual reconstruction network for image compressive sensing,” Neurocomputing 359, 483–493 (2019).
  * [41] K. Kulkarni, S. Lohit, P. Turaga, R. Kerviche, and A. Ashok, “Reconnet: Non-iterative reconstruction of images from compressively sensed measurements,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition,_ (2016), pp. 449–458.
  * [42] J. M. Bioucas-Dias and M. A. Figueiredo, “A new twist: Two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Transactions on Image processing 16, 2992–3004 (2007).
  * [43] X. Yuan, Y. Liu, J. Suo, and Q. Dai, “Plug-and-play algorithms for large-scale snapshot compressive imaging,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition,_ (2020), pp. 1447–1457.
  * [44] S. H. Chan, X. Wang, and O. A. Elgendy, “Plug-and-play admm for image restoration: Fixed-point convergence and applications,” IEEE Transactions on Computational Imaging 3, 84–98 (2016).
  * [45] J. Rick Chang, C.-L. Li, B. Poczos, B. Vijaya Kumar, and A. C. Sankaranarayanan, “One network to solve them all–solving linear inverse problems using deep projection models,” in _Proceedings of the IEEE International Conference on Computer Vision,_ (2017), pp. 5888–5897.
  * [46] C. Metzler, A. Mousavi, and R. Baraniuk, “Learned d-amp: Principled neural network based compressive image recovery,” in _Advances in Neural Information Processing Systems,_ (2017), pp. 1772–1783.
  * [47] J. Zhang and B. Ghanem, “Ista-net: Interpretable optimization-inspired deep network for image compressive sensing,” in _Proceedings of the IEEE conference on computer vision and pattern recognition,_ (2018), pp. 1828–1837.
  * [48] J. Sun, H. Li, Z. Xu _et al._ , “Deep admm-net for compressive sensing mri,” in _Advances in neural information processing systems,_ (2016), pp. 10–18.
  * [49] A. Bora, A. Jalal, E. Price, and A. G. Dimakis, “Compressed sensing using generative models,” arXiv preprint arXiv:1703.03208 (2017).
  * [50] Y. Wu, M. Rosca, and T. Lillicrap, “Deep compressed sensing,” arXiv preprint arXiv:1905.06723 (2019).
  * [51] D. Ulyanov, A. Vedaldi, and V. Lempitsky, “Deep image prior,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition,_ (2018), pp. 9446–9454.
  * [52] Y. Wang, L. Lin, Q. Zhao, T. Yue, D. Meng, and Y. Leung, “Compressive sensing of hyperspectral images via joint tensor tucker decomposition and weighted total variation regularization,” IEEE Geoscience and Remote Sensing Letters 14, 2457–2461 (2017).
  * [53] K. M. León-López and H. A. Fuentes, “Online tensor sparsifying transform based on temporal superpixels from compressive spectral video measurements,” IEEE Transactions on Image Processing 29, 5953–5963 (2020).
  * [54] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE transactions on image processing 13, 600–612 (2004).
  * [55] M. Marquez, H. Rueda-Chacon, and H. Arguello, “Compressive spectral light field image reconstruction via online tensor representation,” IEEE Transactions on Image Processing 29, 3558–3568 (2020).
  * [56] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition,_ (2016), pp. 770–778.
  * [57] J. Masci, U. Meier, D. Cireşan, and J. Schmidhuber, “Stacked convolutional auto-encoders for hierarchical feature extraction,” in _International conference on artificial neural networks,_ (Springer, 2011), pp. 52–59.
  * [58] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention,_ (Springer, 2015), pp. 234–241.
  * [59] B. Arad and O. Ben-Shahar, “Sparse recovery of hyperspectral signal from natural rgb images,” in _European Conference on Computer Vision,_ (Springer, 2016), pp. 19–34.
  * [60] A. Chakrabarti and T. Zickler, “Statistics of real-world hyperspectral images,” in _CVPR 2011,_ (IEEE, 2011), pp. 193–200.
  * [61] C. V. Correa, H. Arguello, and G. R. Arce, “Spatiotemporal blue noise coded aperture design for multi-shot compressive spectral imaging,” JOSA A 33, 2312–2322 (2016).
  * [62] L. Galvis, E. Mojica, H. Arguello, and G. R. Arce, “Shifting colored coded aperture design for spectral imaging,” Applied optics 58, B28–B38 (2019).