Open Access
Add to BibSonomy
Abstract
A novel reconstruction method for compressive spectral imaging is designed by assuming that the spectral image of interest is sufficiently smooth on a collection of graphs. Since the graphs are not known in advance, we propose to infer them from a panchromatic image using a state-of-the-art graph learning method. Our approach leads to solutions with closed-form that can be found efficiently by solving multiple sparse systems of linear equations in parallel. Extensive simulations and an experimental demonstration show the merits of our method in comparison with traditional methods based on sparsity and total variation and more recent methods based on low-rank minimization and deep-based plug-and-play priors. Our approach may be instrumental in designing efficient methods based on deep neural networks and covariance estimation.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Compressive spectral imaging (CSI) systems aim at capturing large volumes of spatio-spectral information of a scene of interest from a set of noisy undersampled observations, providing the means to collect such information in settings such as autonomous navigation and exploration, where there may be restrictions on acquisition time and power consumption. In doing so, a noisy system of linear equations (representing the relationship between the spectral image of the scene and its coded snapshot) needs to be solved for a spectral image estimate. Since such a system of linear equations is underdetermined, a fundamental element to solving the problem is to obtain a solution that satisfies a set of pre-specified regularity constraints on the set of spectral images that comply with the measurements.
In the last decade, various approaches to regularize the CSI problem have been proposed, from the most traditional priors based on sparse modeling to the most recent based on deep learning [1,2]. When we assume that a spectral image has a sparse representation on a given dictionary, the solution to the CSI problem is the spectral image with the sparsest representation that fits the measurements to a certain level. Since the sensing matrix often has a structure with limited levels of randomness, there may be a lack of incoherence between the dictionary and the sensing matrix, resulting in loss of reconstruction accuracy. Motived by this premise, various methods based on structured sparsity, non-local means regularization, and low-rank constraints have been proposed, e.g., [3–5]. In [5], for instance, sets of spatio-spectral patches indexed by spatially proximate coordinates are assumed to form submatrices with low-rank structure, and therefore by denoising such submatrices via weighted nuclear norm minimization in an iterative fashion, it is possible to obtain highly accurate spectral image estimates.
In addition to acquisition speed and energy-efficient transmission of spectral images, new emerging applications also demand expeditious reconstructions [2,6]. Due to the availability of large datasets and the advances in deep-based architectures, recently there is a lot of interest in developing fast and efficient learning-based reconstruction methods to tackle inverse problems in general [2,7,8]. This has led to the development of end-to-end convolutional networks for spectral image reconstruction using self-attention [9], and most recently deep plug-and-play priors for regularization of CSI problems [6]. From an statistical point of view, these methods exploit, in one way or another, an estimate of the prior distribution $p(x)$ over the set of spectral images, which is inferred from spectral image examples [10,11]. However, this requires a considerable amount of examples, or high sample complexity, particularly in large-scale high dimensional problems such as CSI. Consequently, even though such methods lead to fast and accurate spectral image estimation, they may be limited by their inability to generalize to new unseen data and the shortage of large and rich training datasets [2,6].
A simple approach to reduce sample complexity requirements utilized by deep or learning-based methods (relying on dictionary learning, Gaussian mixture models (GMMs), and end-to-end learning) in the context of CSI is to operate on spatio-spectral patches or blocks of the spectral image as opposed to the entire image [3,12,13]. In this setting, since we can search for block estimates independently, it is possible to attain reconstruction speedups by parallel computing. Here, the solution to the problem often results from merging a collection of local estimates, and as a result such speedups may come at the expense of global optimality. In practice, however, this may not be too problematic because we often have access to noisy measurements and an imperfect imaging model, thereby limiting the overall global performance in any case.
The idea of this work is to exploit a block-based reconstruction framework, but from a different perspective. Inspired by recent advances in graph signal processing (GSP), we rely on the notion of smoothness on a graph to design a simple yet efficient signal reconstruction approach for CSI. Specifically, when the signal of interest is assumed to be bandlimited on a graph, an efficient signal estimate with a desirable closed-form can be obtained by minimizing a differentiable function, referred to as the graph Laplacian quadratic form provided a few samples of the signal on the graph [14,15]. In our setting, an important question to ask is: what graph should we use to accurately reconstruct spectral images? the answer is that such graph must be inferred from a suitable set of signals residing on the graph in such a way that the subject signals are the smoothest on the graph. Ultimately, the hope is that the original signal is sufficiently smooth on such a graph so that an accurate signal estimate can be obtained. The contributions of this work are as follows.
- 1. We show that given a signal obtained from a panchromatic image of the scene, a suitable graph can be inferred, leading to accurate spectral image estimates.
- 2. We propose a method to infer the graph, which is based on a recent state-of-the-art graph learning method, and also we design a block-based reconstruction approach based on the inferred graphs, leading to notable performance against traditional and modern methods.
- 3. We provide an extensive simulation study and an experimental demonstration of our approach using a single disperser coded apertature snapshot spectral imaging (SD-CASSI) system with a dual camera.
- 4. We discuss aspects such as graph selection and spatial spectral resolution, provide theoretical principles to define the block size as well as limitations of our current approach.
We note that unlike dictionary learning, GMMs, and end-to-end learning, which are often learned by means of solving complex non-convex problems provided a significant amount of training examples [16–18], our graph-based method only depends on inferring (or learning) a graph from a panchromatic image using an efficient convex optimization program. This in turn may lead to the creation of fast graph learning methods with desirable optimality guarantees. Moreover, given a suitable graph, a smooth solution on the graph can be found by solving a sparse system of linear equations, which is significantly easier than searching for sparse solutions, or computing an estimate based on the posterior distribution associated with a Gaussian mixture prior model. By designing linear equation solvers specialized to the CSI task, our method arguably can be accelerated as much as end-to-end deep-based methods. In our context, convergence, stability, and accuracy analyses can be performed almost effortlessly with existing mathematical tools, thus allowing a thorough characterization of graph-based methods for CSI. Also, our methods eventually can be applied to CSI systems with no additional cameras. Examples include a dual-disperser CASSI [19,20] with a digital micromirror device, where the side information is attained in a snapshot when all mirrors are on, and the spatial-spectral encoded compressive imager [21,22], where its coded snapshots with little or no spatial blurring may serve as side information. These topics, however important, are out of the scope of the paper.
The remainder of the paper is organized as follows. Section 2. sets the ground for CSI reconstruction using a block-based reconstruction framework. Section 3. describes fundamental graph signal processing concepts, explains how to use the concept of smoothness on graphs in the context of the block-based reconstruction framework, and concludes with the proposed graph-based reconstruction method. Section 4. evaluates the proposed method using simulated and real measurements and compares its performance with traditional and modern approaches. Section 5. discusses important aspects that influence the accuracy of the proposed approach, and Section 6. concludes the paper.
2. Block-based spectral image estimation from compressive measurements
The goal of this section is thus to set the grounds for the design of block-based reconstruction algorithms that rely on graphs for regularization.
Consider a spectral image $X_1,{\ldots},X_L \in \mathbb {R}^{n_1 \times n_2}$ with $L$ bands of size $n_1 \times n_2$, and let $x \in \mathbb {R}^{n_1n_2L}$ such that $x=(\textrm{vec}(X_1)^T,{\ldots},\textrm{vec}(X_L)^T)^T$ be its vector representation. The idea is to reconstruct $x$ from a compressive snapshot $y \in \mathbb {R}^{n_1(n_2+L-1)}$, obeying
where $\mathbf {H}$ is assumed to be a sensing matrix obtained by discretizating the continous imaging model of a single-disperser coded aperture snapshot spectral imager (SD-CASSI) [23], and $e$ is an additive noise term.Recall that the spectral image $x \in \mathbb {R}^{n_1n_2L}$ resides (originally) on a three-dimensional domain given by
that is, the value of $x$ at the entry $k$ is associated with a three dimensional coordinate $(k_1,k_2,k_3) \in \Omega$. Similarly, a CASSI snapshot $y \in \mathbb {R}^{n_1(n_2+L-1)}$ resides on a two-dimensional domain given byDue to the structure of $\mathbf {H}$, for a given measurement subdomain (or subset) $\tilde {\Omega }_k \subset \tilde {\Omega }$ such that
Provided $\tilde {\Omega }_k$, there are multiple $\Omega _k$’s that will satisfy Eq. (2). We thus need to select an appropriate $\Omega _k$. To illustrate, consider a CASSI system with $n_1=n_2=15$ and $L=6$. Given $\tilde {\Omega }_k$ with parameters $b_1=b_2=4$, depicted by the purple circles in Fig. 1, the pair of subsets depicted by the red circles in Fig. 1 both define a suitable $\Omega _k$, and thus satisfy Eq. (2). For the purpose of regularization however, it is best to select $\Omega _k$, taking form of the set in Fig. 1 (left), which can be found as follows:
- 1. Construct a CASSI matrix $\tilde {\mathbf {H}}$ with an all ones coded aperture, i.e., $100$ percent transmittance, as in the Appendix.
- 2. Given $\tilde {\Omega }_k$, solve where $\tilde {\mathbf {H}}_{\tilde {\Omega }_k,\Omega _k^c}$ is the submatrix of $\tilde {\mathbf {H}}$ with rows and columns indexed by $\tilde {\Omega }_k$ and $\Omega _k^c=\Omega \setminus \Omega _k$ respectively. The subset $\Omega _k$ in Fig. 1 (right) results from setting $\tilde {\textbf {H}}$ as a CASSI matrix with 50-percent transmittance binary random coded aperture.
Fig. 1. Illustration of pairs of subsets $(\tilde {\Omega }_k,\Omega _k)$ satisfying Eq. (2). Here $\tilde {\Omega }_k$ is fixed, and $\Omega _k$ can be selected as either the set of red circles on the left or on the right. In this work, we define pairs $(\tilde {\Omega }_k,\Omega _k)$ using Eq. (3), where $\Omega _k$ takes the form depicted on the left.
We now formally state an algorithm for block-based spectral image estimation. Define a collection of overlapping subsets $\tilde {\Omega }_k \subset \tilde {\Omega }$ as in Eq. (1) such that
where $\mathcal {K}$ denotes a subset of $\tilde {\Omega }$ that makes adjacent $\tilde {\Omega }_k$’s be overlapped to some degree. In this work, $\mathcal {K} \subset \tilde {\Omega }$ is given by:Next, find the collection $\Omega _k \subset \Omega$ associated with $\tilde {\Omega }_k$ by using Eq. (3). Given a collection of convex regularity functions $\mathcal {R}_k : \Omega _k \mapsto [0,\infty )$, a block-based estimate $\hat {x}$ of a spectral image can be obtained as follows:
- 1. For each pair $(\tilde {\Omega }_k,\Omega _k)$, solve$$\hat{x}_k = \arg\min_{x \in \mathbb{R}^{|\Omega_k|}} \ \mathcal{R}_k(x) \ \text{s.t.} \ \|\mathbf{H}_{\tilde{\Omega}_k,\Omega_k}x - y_{\tilde{\Omega}_k}\|_2^2\leq \epsilon_k^2,$$where $\epsilon _k$ is such that the noise term $e_k$ satisfies $\|e_k\|_2\leq \epsilon _k$.
- 2. Given the local estimates $\hat {x}_k$, the estimates can be combined into a global estimate $\hat {x}$ by solving the problem where $\mathbf {S}_k \in \{0,1\}^{|\Omega _k|\times n_1n_2L}$ is such that $\mathbf {S}_k x$ outputs the subvector $x_{\Omega _k}$ of $x$ indexed by $\Omega _k$ for all $x \in \mathbb {R}^{|\Omega |}$.
We note that Eq. (6) is a common strategy to combine local estimates in several image processing tasks such as image restoration based on group-sparse representations [25], and image fusion based on Gaussian mixtures [26] to name a few. It is easy to show that the solution to Eq. (6) has closed-form given by
where the matrix $\sum _{k}\mathbf {S}_k^T\mathbf {S}_k$ can proven diagonal, and its inverse can therefore be computed efficiently.3. Smoothness on graphs for block-based spectral image estimation
3.1 GSP preliminaries
An undirected graph $G=(V,E,w)$ is a triple, consisting of a vertex set $V=\{v_i\}_{i=1}^n$, an edge set $E \subset V \times V$, and a non-negative weight function $w: E \mapsto [0,\infty )$ s.t. $w(i,j)=w(j,i)$ for $(i,j) \in E$, and $w(i,j)=0$ for $(i,j) \notin E$. Also, it is assumed that the graph has no self-loops and thus $(i,i) \notin E$ for $i \in V$.
A graph signal $x$ on $G$ is a function $x: V \mapsto \mathbb {R}$ such that the value of $x$ at $i \in V$ is given by $x_i \in \mathbb {R}$.
The adjacency matrix $\mathbf {W}_G$ of $G$ is an $n \times n$ matrix whose entries $W_{ij}$ are given by the weight function $w(i,j)$ at the edge $(i,j) \in V \times V$, i.e., $W_{ij}=w(i,j)$. By definition $w(i,j)=w(j,i)$, and therefore we have that the adjacency matrix $\mathbf {W}_G$ is symmetric, i.e., $\mathbf {W}_G=\mathbf {W}^T_G$. The degree matrix $\mathbf {D}_G$ of $G$ is an $n \times n$ diagonal matrix whose diagonal entries $D_{ii}$ are given by $\sum _{j=1}^n w(i,j)$. The graph Laplacian $\mathbf {L}_G \in \mathbb {R}^{n \times n}$ of $G$ is an $n \times n$ matrix whose diagonal entries $L_{ii}$ are given by $\sum _{j=1}^n w(i,j)$ and off-diagonal entries $L_{ij}$ are given by $-w(i,j)$, equivalently
An important characteristic of the graph Laplacian is that it is symmetric positive semidefinite, thus its eigenvalues are real and nonnegative.3.2 Smoothness on graphs and learning graphs from data
Under the assumption that the spectral image $x$ at $\Omega _k$ is smooth with respect to a graph $G_k$ with vertex set $\Omega _k$, the regularization function $\mathcal {R}_k(x)$ in Eq. (5) can be defined as:
where $\mathbf {L}_{G_k}$ is the graph Laplacian of $G_k$.Since the $G_k$’s are not known a priori, we use a co-registered panchromatic image of the scene $Z \in \mathbb {R}^{n_1 \times n_2}$ to infer the structure of the graphs. From a GSP perspective, the problem of learning graphs from data requires a collection of signals residing on the vertices of the graph to be inferred and there are numerous methods to do so as recently surveyed by [27,28]. Here, we adapt an efficient convex optimization program to the CSI context, recently proposed by [29], that has nice scalability properties on the number of vertices of the graph. To begin, we partition the signal sub-domain $\Omega _k = \bigcup _{l=l_0}^{L'}\Omega _k^l$ with $1 \leq l_0 \leq L' \leq L$ in such a way that
and define a signal $z_k \in \mathbb {R}^{|\Omega _k|}$ on $\Omega _k$ based on the panchromatic image $Z \in \mathbb {R}^{n_1 \times n_2}$ as follows:The goal is thus to construct a weighted graph $G_k$ with vertex set $\Omega _k$ and adjacency matrix $\mathbf {W}_{G_k}$ by using the statistical properties of $z_k$, i.e., radiometric distances among the elements of $z_k$. To do so, let $\mathcal {W}_{|\Omega _k|}$ be a set of admissible adjacency matrices such that
The former edge constraint says that there is no relation between pixels that belong to different spectral bands, and therefore it seeks to impose spatial regularity on each band separately. The latter edge constraint on the other hand limits the family of graphs to that of q-nearest neighbor graphs, and controls the sparsity of the resulting graph. The solution to Eq. (8) can be found by using the publicly available algorithm proposed by Kalofolias et al. [29,30].
Algorithm 1 summarizes the block-based spectral image estimation procedure based on graph constructions resulting from Eq. (8).
Algorithm 1. Our method
We note that Step 6 of Algorithm 1 is equivalent to solving $\textrm{min}_{x_k\in \mathbb {R}^{|\Omega _k|}} \ x_k^T\mathbf {L}_{G_k} x_k \ \textrm{s.t.} \ \mathbf {H}_{\tilde {\Omega }_k,\Omega _k}x_k=y_{\tilde {\Omega }_k}$ when $\epsilon _k=0$, which has a solution $\hat {x}_k$ obeying:
Also, the selection of $\theta$ and $q$ is fundamental to finding suitable graphs that encode structural information from $z_k$. For simplicity, in our experiments, we assume $\theta$ and $q$ to be the same across different patches. However, it may well be that we need to adjust them according to the patch characteristics for improved results. Similarly, according to our preliminary experiments, we notice that the value of $\theta$ depends on $q$. In this work, we did not pay special attention to the fine-tuning of such parameters, but these are nonetheless important for optimal performance.
4. Experimental results
We now evaluate the reconstruction performance of the proposed Algorithm 1 to reconstruct spectral images from simulated and real measurements. In particular, we assume the CSI system to be single disperser coded aperture snapshot spectral imaging (SD-CASSI) system. The performance of our approach is compared with traditional reconstruction approaches using total variation TV and sparsity WDCT and more recent methods based on low-rank minimization DeSCI [5] and deep-based plug-and-play priors PnP [6]. TV searches for a spectral image estimate by minimizing the total-variation norm $\|x\|_{\textrm{TV}}=\sum _{l=1}^L\sum _{i,j}\|\nabla (X_l)_{ij}\|_2$ over the set of spectral images defined by the measurements, i.e., $\{x:\|\mathbf {H}x-y\|_2^2\leq \epsilon ^2\}$, where $\epsilon$ denotes the noise level, and $\nabla \approx (\frac {d}{dx},\frac {d}{dy})^{\textrm{T}}$ is a first order finite difference approximation of the gradient [31]. WDCT searches for a spectral image estimate by minimizing $\|\mathbf {\Psi }^Tx\|_1$ over the set of signals defined by the measurements, where $\mathbf {\Psi } = \mathbf {\Psi }_{\textrm{2D-W}} \otimes \mathbf {\Psi }_{\textrm{1D-DCT}}$ such that $\mathbf {\Psi }_{\textrm{2D-W}}$ denotes two-dimensional Symmlet-8 wavelet transform basis, and $\mathbf {\Psi }_{\textrm{1D-DCT}}$ denotes discrete cosine transform basis [1,32–34].
Since our approach relies on side information to construct the graph, in simulations we adapted the TV and WDCT approaches to consider side information as a direct measurement. This means that the both the CASSI’s sensing matrix $\mathbf {H}$ and the side information sensing matrix $\mathbf {R}$ define the feasible set of signals, i.e., $\{x \in \mathbb {R}^{n_1n_2L}: \mathbf {H}x=y,\mathbf {R}x=z\}$ in the noiseless case, $\epsilon =0$. In practice, however, this requires an extra calibration step of estimating the side information matrix, and thus in the case of real measurements, we don’t use the side information as a direct measurement. This also illustrates that our approach can exploit the side information even if we do not know such a matrix in advance.
4.1 Simulated measurements
We consider spectral images $\{X_l\}_{l=1}^L$, coming from the spectral image databases: UDEL in Fig. 2, HRVRD [35], and ICVL [36]. Particularly, HRVRD and ICVL refers to a collection of eight selected datacubes with diverse spatial spectral content, where HRVRD consists of img1, imga6, imgb3, imgb8, img3, imga8, imgd6, imgg8; and ICVL consists of gavyam_0823-0930, gavyam_0823-0933, gavyam_0823-0944, plt_0411-1037, plt_0411-1211, plt_0411-1232-1, rsh_0406-1356, rsh_0406-1413. For simplicity, we downsampled all spectral images to size $n_1 \times n_2 \times L$ with $n_1=n_2=512$ and $L=31$ and also normalized by their maximum intensity across the spatial-spectral coordinates. The goal is to reconstruct $\{X_l\}_{l=1}^L$ from a single snapshot $Y \in \mathbb {R}^{n_1 \times n_2 + L -1}$ satisfying:
where again $x=(\textrm{vec}(X_1)^T,{\ldots},\textrm{vec}(X_L)^T)^T$, $y=\textrm{vec}(Y)$, and $\mathbf {H}$ denotes a single disperser CASSI matrix with binary random coded aperture of size $n_1 \times n_2$ and $50 \%$ transmittance, constructed as in the Appendix .Fig. 2. UDEL hyperspectral image database collected at our lab. RGB renderings of four hyperspectral datacubes (HSDC) used for performance evaluation of the signal recovery algorithms. The hyperspectral datacubes HSDC1, HSDC2, HSDC3, and HSDC4 consist of 31 spectral bands of size 2064-by-3088 pixels, ranging from 400 to 700 nm.
Provided a panchromatic image $Z \in \mathbb {R}^{n_1 \times n_2}$ such that $Z=\frac {1}{L}\sum _{l=1}^LX_l$, we estimate $x$ from $y$ using the proposed block-based approach, described in Sec. 3. with sub-domain parameters $b_1=b_2=32$ and $\mathcal {K}$ as in Eq. (4) and graph-learning parameters $\theta =100$, $q=33$. As aforementioned, TV and WDCT were run in two settings with w$/$ and without w$/$o side information as a direct measurement. Spectral image reconstructions with TV and WDCT were obtained using C-SALSA [37], and the noise level $\epsilon$, the maximum number of iterations MAXIT, and the tolerance TOL were set to $10^{-8}$, $10000$, and $10^{-6}$, respectively. Recall that the noise level $\epsilon$ determines the extent of the feasible set, so by setting it with a small enough value, it is suggested that the measurements are noiseless. In the case of DeSCI and PnP, reconstructions were obtained by the publically available code provided by the authors of [5,6] with the default settings. Lastly, our spectral image reconstructions were obtained by using Algorithm 1, where Step 6 was solved using the MINRES algorithm (proposed in [38]) with the tolerance TOL and the maximum number of iterations MAXIT equal to $10^{-6}$ and $10000$ respectively.
Table 1 shows the average reconstruction performance across the selected databases. Observe that the proposed method compares notably against all methods, however TV w$/$ performs slightly better than our approach in terms of mean structural similarity index (SSIM), mean peak signal-to-noise ratio (PSNR), and mean spectral angle mapper (SAM) as defined in [39,40]. This can be explained due to the fact that our block-based method does not exploit the global structure of compressive measurements to narrow down the space of spectral images, as opposed to TV w$/$ that searches for image estimates in the set $\{x \in \mathbb {R}^{n_1n_2L}: \mathbf {H}x=y,\mathbf {R}x=z\}$. Our approach may thus lose accuracy in simulations. In practice, however, this may not be too problematic because we rarely want to fit the measurements to machine precision; not only is the sensing matrix imperfect but the measurements are also noisy. Without knowing the side information matrix $\mathbf {R}$ in advance, our approach is capable of attaining highly accurate spectral image estimates. To illustrate, Fig. 3 displays reconstructed RGB renderings of the HSDC1, obtained using our approach and the rest of the approaches without side information as a direct measurement.
Fig. 3. Reconstructions from a simulated CASSI snapshot. RGB rendering of the original spectral datacube, and reconstructed RGB renderings obtained by TV w$/$o, WDCT w$/$o, DeSCI, PnP, and Our method with $b_1=b_2=32$, where w$/$o indicates that side information was not used as a direct measurement. We should note, however, that our graph-based method exploits statistical information from a panchromatic image to infer the graphs. Also in the figure, original and reconstructed spectral signatures at locations P1, P2, and P3 are displayed on the original RGB image. In the legends, the values in parenthesis are SAM values in degrees; the smaller the better.
Table 1. Average performance metrics of our method and comparing approaches TV, WDCT, DeSCI, PnP across different databases UDEL(4), HRVRD(8), and ICVL(8), where the parenthesis encloses the number of datacubes per database. w$/$ and w$/$o indicate whether the reconstruction algorithm uses side information or not as a direct measurement of the spectral image of interest. We note however that our method still uses the side information image to construct the graphs. Our method’s R-TIME includes in parenthesis average overhead times due to graph construction and aggregation, GC-TIME and A-TIME.
Also in Table 1, we show the reconstruction time R-TIME of all algorithms except the PnP’s since we didn’t obtain an accurate estimate of the PnP’s R-TIME. The reason is that PnP reconstructions were obtained on a different computing platform while executing several other workloads, leading to longer reconstruction times than those expected in the literature. Note that we did not optimize the maximum number of iterations for the proposed and traditional methods. Instead, these iterated until their solution converged within a given tolerance TOL. We can thus speed up reconstruction in all cases by tunning the number of iterations. We can observe that, in a best-case scenario where each block is processed in parallel, our approach is several orders of magnitude faster than the comparing methods even though there is a little overhead due to graph learning time GC-TIME and aggregation time A-TIME, which, we expect, can be substantially decreased by code optimization and more efficient graph learning. To clarify, for our approach, the reported R-TIME is the average reconstruction time of a computing system with as many cores as blocks so each block can be processed in parallel. This gives a relation between blocks and cores $C=\#\textrm{Blocks}/\textrm{\#Cores}=1$. In this setting, we would need to have $\textrm{\#Cores}=|\mathcal {K}|=1088$ in order to get $C=1$.
4.2 Real measurements
We now illustrate the performance of our approach in a real experimental setting. In doing so, we implemented the dual-camera compressive spectral imaging (DC-CSI) system in Fig. 4. Our DC-CSI system consists of an objective lens, CoastalOpt UV-VIS-IR 60 mm Macro, that images the scene of interest onto a bandpass filter, which limits the spectral range of the scene to 450-650 nm. The filtered image is then relayed onto the CASSI’s optical path and the side-information camera by a 4-F system with a beam splitter located at the Fourier plane.
Fig. 4. Implemented dual-camera compressive spectral imager. Top right corner displays an optical diagram of the system. The interested reader is referred to [41] for more details about the similar system’s working principle.
In the CASSI’s arm, the image is spatially coded by a random coded aperture of size $128 \times 128$ with $50 \ \%$ transmittance, and next the spatially-coded image is dispersed horizontally by a double Amici prism. Lastly, the spatio-spectrally coded image is captured by the focal plane array FPA1; a camera Basler ac A3088-57um with resolution $3088 \times 2064$ with pixel size $2.4 \times 2.4 \ \mu m$. We set the relay lens such that the $128 \times 128$ coded aperture pixels were mapped onto about $1045 \times 1045$ camera pixels. Since the coded aperture pixels are about eight times as big as the FPA1 pixels, the size of the coded image is almost unaffected by the relay lens.
In the side-information arm of the system, a panchromatic image of the scene is captured by the focal plane array FPA2; a camera Basler ac A3088-57uc with resolution $3088 \times 2064$ and pixel size $2.4 \times 2.4 \ \mu m$. Such an image was optically aligned to the $550~nm$ band captured by the FPA1. Due to slight variations in the image sizes captured by both the FPA1 and FPA2, we noticed some degree of alignment error (2-3 pixels), which may lead to some loss of reconstruction accuracy. We assume, however, perfect alignment when performing reconstructions. Table 2 briefly summarizes the main specifications of our system.
Table 2. System specifications of our DC-CSI system, displayed in Fig. 4.
The DC-CSI system was configured to encode $26$ spectral bands of size $256 \times 256$ from a scene of interest onto a compressive snapshot of size $256 \times 281$. In doing so, given a desired spatial resolution, the number of spectral bands and associated wavelengths were determined by an estimation of the prism dispersion curve. Then, provided the selected wavelengths, the response of the system at each wavelength (often referred to as the calibration datacube) was obtained. The interested reader is referred to [23,42] for more details about this process.
Given the compressive and panchromatic snapshots, we assume that the spectral image of interest obeys:
where $\tilde {y}:=\textrm{vec}(Y)/\max _{i,j}(Y_{ij})$, $\tilde {e}$ is an additive noise term, and $\hat {\mathbf {H}}$ represents the sensing matrix of our DC-CSI system, constructed from the calibration datacube as described in the Appendix. Spectral image estimates of SCN1 and SCN2 in Fig. 5 were then obtained by solving 11 using the approaches TV w$/$o, DeSCI, PnP, and our method with sub-domain parameters $b_1=b_2=32$ and $\mathcal {K}$ as in Eq. (4), and graph-learning parameters $\theta =100$, $q=33$.Fig. 5. Color images of the spectral scenes of interest. The color images were captured by the side-information camera in Fig. 4 under the same illumination conditions used for the CASSI snapshot. The spectral signatures at the points P1, P2, and P3 will be evaluated.
An important parameter to enable block-based reconstructions with real data is the noise level $\epsilon _k$. Since the problem Eq. (5) can be reformulated as:
To simplify notation, let $\tilde {y}_k$ and $\hat {\mathbf {H}}_k$ represent $\tilde {y}_{\tilde {\Omega }_k}$ and $\hat {\mathbf {H}}_{\tilde {\Omega }_{k},\Omega _{k}}$ respectively. Furthermore, let $k^*$ index the patch $\tilde {y}_{k^*}$ with the highest variance over $\mathcal {K}$, that is $k^*=\arg \max _{_{k\in \mathcal {K}}} \ \textrm{SD}(\tilde {y}_{k})$, where $\textrm{SD}(\tilde {y}_{k})$ computes the standard deviation of the entries of $\tilde {y}_k$.
- 1. For each $\alpha \in \textrm{logspace}(-2,2,50)$, solve for $x_{\alpha }$ using Eq. (12) with $k:=k^{*}$ and $\alpha _k:=\alpha$, and record the values ($\|\tilde {y}_{k^*} - \hat {\mathbf {H}}_{k^*}x_{\alpha }\|_2^2$, $x^T_{\alpha }\mathbf {L}_{G_{k^*}}x_{\alpha }$)
- 2. Set $\alpha ^*$ as$$\arg\min_{\alpha} \ |C_1^{{-}1}\|\tilde{y}_{k^*} - \hat{\mathbf{H}}_{k^*}x_{\alpha}\|_2^2-C_2^{{-}1}x^T_{\alpha}\mathbf{L}_{G_{k^*}}x_{\alpha}|,$$where $C_1=\max _{\alpha }\|\tilde {y}_{k^*}-\hat {\mathbf {H}}_{k^*}x_{\alpha }\|_2^2$ and $C_2=\max _{\alpha }x_{\alpha }^T\mathbf {L}_{G_{k^*}}x_{\alpha }$.
We used the just explained procedure to fine-tune the noise level parameter $\epsilon$ for TV, and the parameters for DeSCI and PnP were adjusted to the best of our capabilities. In PnP, it was crucial to set correctly the noise level sequence, which has to decrease as the iteration counter progresses, and the maximum number of iterations.
The accuracy of spectral image estimates in the visible spectrum can be evaluated to some extent by looking at their RGB representation. Figures 8 and 9 show reconstructed RGB images of SCN1 and SCN2 obtained through TV w$/$o, DeSCI, PnP, and our method from the measurements in Fig. 6. Here, TV w$/$o, DeSCI, PnP operate under their nominal conditions. That is, they do not exploit side information as a direct measurement while our approach uses the side information only to construct the required collection of graphs. We can observe in the figure that our method produces estimates with as many spatial details as those in the side information and also with as many colors as those present in the original color image of the scenes. The reconstructed RGB images exhibit, however, chromatic aberrations due to imperfect separation of the spectral bands from the coded snapshot. This results in color smearing, clearly seen in the TV-based reconstruction, where colors from nearby objects mix, forming new colors at edge transitions. Specifically, our approach may lead to solutions with block artifacts when the regularized sensing matrix, $\mathbf {H}^T_k\mathbf {H}_k + \alpha \mathbf {L}_{G_k}$ is non-invertible or singular. This is due to the nature of the null space of the graph Laplacian matrix, which does not penalize certain piece-wise constant signals, and they thus may appear in the solution without affecting the overall cost of the objective function.
Fig. 6. Compressive and panchromatic snapshots of the scenes SCN1, SCN2 as seen by our DCCSI system displayed in Fig. 4.
Fig. 7. From left to right, selection of regularization parameter $\alpha ^*$ for SCN1 and SCN2 in Fig. 5, respectively. The patch size is determined by $b_1=b_2=32$.
Fig. 8. Reconstructions from real measurements. RGB image of the original scene obtained by a color camera, and reconstructed RGB renderings obtained by our method, TV, DeSCI, and PnP without side information as direct measurement. Also in the figure, reconstructed spectral signatures at a few locations P1, P2, and P3 displayed in ORIG. In the legends, the values in parenthesis indicate the accuracy of the reconstruction in terms of the Spearman’s correlation coefficient, which is a robust measure of association between two variables.
Fig. 9. Reconstructions from real measurements. RGB image of the original scene obtained by a color camera, and reconstructed RGB renderings obtained by our method, TV, DeSCI, and PnP without side information as direct measurement. Also in the figure, reconstructed spectral signatures at a few locations P1, P2, and P3 displayed in ORIG. In the legends, the values in parenthesis indicate the accuracy of the reconstruction in terms of the Spearman’s correlation coefficient, which is a robust measure of association between two variables.
Also in Figs. 8 and 9, we compare the reconstructed spectral signatures at a few points on the scenes with reference spectra, measured by a non-imaging spectrometer. The reconstructed spectra were generated by averaging the intensities in a small window around a certain point for each band. Both estimated and reference spectral signatures were then normalized by the sum of their intensities. In the figure, we can observe that the reconstructed spectral signatures are strongly correlated to the reference spectra, however, all methods produce outlying variations to some extent at the ends of the spectrum. Our approach and PnP lead to more spectrally accurate reconstructions than the rest of the methods. We should note nonetheless that the PnP seems to handle better the outlying variations at the ends of the spectrum, but our approach seems to outperform at detecting major peaks and transitions along the spectrum.
Figures 10 and 11 allow one to get a closer look into how the reconstructed spectral bands vary for each reconstruction approach. We can observe TV produces spectral images with subtle differences in the spatial content across the spectrum while our approach, DeSCI and, PnP appear to have more distinguishable changes across the spectrum. We can also see that our approach and PnP produce similar results however our approach can preserve better spatial details because our method uses a scene adapted graph-based prior as opposed to PnP which uses a class adapted prior, that is a prior trained on the class of spectral images.
Fig. 10. Spectral image estimates of the scene SCN1 in Fig. 5. We display 8 out of 26 reconstructed spectral bands of size $256 \times 256$. From left to right, the TV, our method, DeSCI, and PnP along the columns.
Fig. 11. Spectral image estimates of the scene SCN2 in Fig. 5. We display 8 out of 26 reconstructed spectral bands of size $256 \times 256$. From left to right, the TV, our method, DeSCI, and PnP along the columns.
5. Discussion
5.1 On the selection of the graph
So far, we have presented a possible way to construct the collection of graphs $G_k$, which is needed to enable our block-based approach for spectral image reconstruction. Since a panchromatic image provides mainly spatial information of the scene, we assume the edge structure of the graph is disconnected across the vertices associated with different spectral bands, expecting that the spatio-spectral information of the scene encoded in the compressive snapshot is sufficient to resolve an accurate estimate of the spectral image. Notably, the results point to accurate spectral image reconstructions, validating our approach for a variety of scenes.
However, we note that when the measurements are poorly specified, a disconnected graph may not properly regularize the block-wise CSI problem. Thus, new graph-based methods must explore or develop graph constructions, which exploit spatial-spectral associations between pixels of different spectral bands. This in turn may lead to more stable and accurate graph-based spectral image estimates. In the literature, there is a whole body of work, e.g., [4,9,43,44], which has successfully used spatial-spectral correlation estimates to enhance CSI reconstructions. It is thus reasonable to suggest that by exploiting similar principles in our context, we can obtain performance gains. We note that a first approach to exploit such a concept in our framework is to impose adjacency matrices with block tridiagonal structure in Eq. (8) (as opposed to block-diagonal structure), and infer them from trichromatic or RGB side information, as an alternative to panchromatic side information. This is because RGB side information contains more information about the spatial-spectral correlation of a target scene than a panchromatic image. The quest of such graph-based models represents an important research opportunity, but it is out of the scope of the paper.
5.2 On the selection of the patch size
Despite the benefits of block-based reconstructions, their performance can be affected by the choice of the patch size. This paper does not analyse the effects of this parameter deeply, but we now give some principled ideas to shed light onto how to select the patch size.
To simplify notation, let $\mathbf {H}_{\Omega _k}$ represent $\mathbf {H}_{\tilde {\Omega }_k,\Omega _k}$. Recall that the fundamental unit of the proposed block-based reconstruction algorithms is given by the system of liner equations in Eq. (10). Here, by definition, the patch size $|\tilde {\Omega }_k|=b_1b_2$ determines the block size $|\Omega _k|$.
We note that for the system of equations in Eq. (10) to have a unique solution, the matrix $\mathbf {H}_{\Omega _k}^T\mathbf {H}_{\Omega _k} + \alpha \mathbf {L}_{G_k}$ must be non-singular, which is the same as saying that the null space of both matrices intersect at the zero vector $\mathcal {N}_{\mathbf {H}_{\Omega _k}} \cap \mathcal {N}_{\mathbf {L}_{G}^{1/2}} = \{\mathbf {0}\}$. Using this condition, the following property regarding the number of measurements can be derived.
Assume $\mathcal {N}_{\mathbf {H}_{\Omega _k}} \cap \mathcal {N}_{\mathbf {L}_{G}^{1/2}} = \{\mathbf {0}\}$, and let $\mathbf {H}_{\Omega _k} \in \mathbb {R}^{|\tilde {\Omega }_k| \times n}$. Then the number of measurements $|\tilde {\Omega }_k|$ must obey
To clarify, note that since $\mathcal {N}_{\mathbf {H}_{\Omega _k}}$ and $\mathcal {N}_{\mathbf {L}_{G}^{1/2}}$ are complementary subspaces, the dimension of their sum $\textrm{dim}(\mathcal {N}_{\mathbf {H}_{\Omega _k}}+\mathcal {N}_{\mathbf {L}_{G}^{1/2}})$ can be written as $\textrm{dim} \ \mathcal {N}_{\mathbf {H}_{\Omega _k}}+\textrm{dim} \ \mathcal {N}_{\mathbf {L}_{G}^{1/2}}$, and is strictly upper bounded by $n$. On the other hand, the theorem of conservation of dimension for matrices [45] gives us that $\textrm{dim} \ \mathcal {N}_{\mathbf {H}_{\Omega _k}} + \textrm{dim} \ \mathcal {R}_{\mathbf {H}_{\Omega _k}} = n$, so since $\textrm{dim} \ \mathcal {R}_{\mathbf {H}_{\Omega _k}} \leq |\tilde {\Omega }_k|$, we get the lower bound $n-|\tilde {\Omega }_k| \leq \textrm{dim} \ \mathcal {N}_{\mathbf {H}_{\Omega _k}}$. We thus obtain that $n > \textrm{dim} \ \mathcal {N}_{\mathbf {H}_{\Omega _k}}+\textrm{dim} \ \mathcal {N}_{\mathbf {L}_{G}^{1/2}}\geq n-|\tilde {\Omega }_k| + \textrm{dim} \ \mathcal {N}_{\mathbf {L}_{G}^{1/2}}$, and Eq. (13) follows.
Since by construction, $G_k$ has at least $L'-l_0+1$ connected components with $1 \leq l_0 \leq L' \leq L$, the dimension of the $\textrm{dim} \ \mathcal {N}_{\mathbf {L}_{G_k}^{1/2}}$ is lower-bounded by $L'-l_0+2$. When $\Omega _k$ forms a parallelepiped as in Fig. 1(left), the patch size $|\tilde {\Omega }_k|$ must satisfy
In practice, the size of $|\tilde {\Omega }_k|$ may grow several times greater than $L+1$ before we can solve the above system in an stable manner. At a qualitative level, we can say that the presence of block artifacts in the solution are an indication that the system is poorly conditioned and therefore unstable. A simple strategy to avoid block artifacts that was put to practice in the experiments with real data and does not require the system to be perfectly stable is to relax the tolerance parameter of the linear equation solvers. This means that we don’t seek to have the exact solution but an approximate one that is good enough. Lastly, recall that stability follows from the fact $\mathcal {N}_{\mathbf {H}_{\Omega _k}} \cap \mathcal {N}_{\mathbf {L}_{G_k}^{1/2}}=\{0\}$. For this condition to hold it is sufficient that the smallest eigenvalue of $\mathbf {H}_{\Omega _k}^T\mathbf {H}_{\Omega _k} + \alpha \mathbf {L}_{G_k}$ is greater than zero. As a result, we can try to characterize such an eigenvalue, and select the patch size in such a way that on average $\tilde {\Omega }_k$ leads to (approximately) well conditioned systems of linear equations.5.3 On the spatial spectral resolution
Unlike the traditional SD-CASSI’s spatial resolution, limited by the coded aperture’s pixel size, our dual camera SD-CASSI system using graphs enables the transfer of the spatial resolution of the side-information camera to the reconstructed spectral image as suggested by our real-data experiments. This is because our system has higher resolution detectors with smaller pixel sizes than the coded aperture’s pixel size. We can thus construct the graph on a higher-resolution image domain, allowing the search of spectral image estimates with higher spatial resolution than the coded aperture’s. On the other hand, in the traditional SD-CASSI system, the spectral resolution is ideally determined by the prism dispersion across the detector’s horizontal direction and the coded aperture’s pixel size [23,42]. In our case, due to the extra spatial resolution enabled by the graph through the side camera, our method could provide some additional spectral resolution as well. We note that we do not intend to give a definite characterization of the spatial-spectral resolution of our system, and a thorough study of it should also consider the effect of the required regularization assumptions on the reconstructed spectral images.
6. Conclusion
We have proposed a novel block-based reconstruction approach to CSI, which is based on the notion of smoothness on graphs, where the graphs are inferred from panchromatic side information. Our approach leads to efficient solutions with a desirable closed-form that can be found by solving multiple sparse systems of linear equations in parallel. By increasing slightly SD-CASSI’s hardware complexity, the new framework can produce notable results with high spatial resolution and minor spectral distortion when compared against traditional and state-of-the-art methods. The proposed graph-aided SD-CASSI system enables both higher spatial and spectral resolutions than traditional SD-CASSI’s as demonstrated by our real data results. When inferring the graphs from a panchromatic image, it is important to be aware that constant regions may lead to graphs that do not represent the image structure suitably because the associated radiometric distances are uninformative. For such cases, it is thus important to design graph learning methods that resort to domain knowledge of the spectral image domain instead of the panchromatic image itself. Since the proposed family of graphs does not exploit spatial-spectral associations among vertices of the spectral image domain (which may lead to more accurate and stable reconstructions), there is room for new graph-based models to arise from our work. Similarly, our results may be instrumental in the design of more efficient deep-based reconstruction methods since it has been shown recently that inductive biases introduced through graph-based smoothness lead to data-efficient and robust deep neural networks for classification tasks [46]. In addition, due to recent advances in covariance matrix estimation from CSI measurements [20], graph-based methods that only rely on the compressive measurements can be designed by casting the graph inference problem on a statistical graph learning framework.
Appendix. Construction of the CASSI sensing matrix
Algorithm 2 explains how to construct the CASSI sensing matrix provided a set of calibration images. The algorithm assumes that the shadow of the coded aperture moves from right to left in steps of one pixel as the CASSI system is (almost) uniformly illuminated with monochromatic light at the $L$ distinct calibration wavelengths. This is suggested by the definition of $\Omega _l$ in line 7 of the algorithm. If the coded aperture’s shadow, on the contrary, moves from left to right in steps of one pixel. Line 7 and 8 have to be updated by replacing $x-L+1$ highlighted in blue with $x-l+1$. Note that the algorithm uses the function $\textrm{sparse()}$ in Line 10, adopted from MATLAB, to construct a sparse matrix of size $n_1(n_2+L-1) \times n_1n_2L$ with nonzero entries determined by $I,J,H$.
For simulations, the calibration images are defined as follows. Consider a random coded aperture $T \in \{0,1\}^{n_1 \times n_2}$ with entries $T_{ij}$ drawn from a Bernoulli distribution with parameter $p=0.5$. The calibration images $\tilde {h}_l$ are defined by
for $(y,x) \in \{1,{\ldots},n_1\} \times \{1,{\ldots},n_2+L-1\}$, where $h_0:\{1,{\ldots} n_1 \}\times \{1,{\ldots} n_2+L-1\} \mapsto [0,\infty )$ such that $h(i,j)=T_{i,j}, \ 1,1 \leq y,x \leq n_1,n_2$, and $h(i,j)=0$ otherwise. One way to make sure that the CASSI matrix has been properly constructed is to verify the following property. Let $\delta _l \in \{0,1\}^{n_1n_2L}$ be uniform monochromatic light at the $l$-th spectral band such that $\delta _l = \textrm{kron}(e_l^T,\textrm{ones}(1,n_1n_2))^T$, where $e_l$ is the $l$-th standard basis vector of size $L \times 1$. Then the CASSI matrix satisfies $\mathbf {H}\delta _l = \textrm{vec}(\tilde {h}_l), \ l=1,{\ldots},L$.Funding
National Science Foundation (1815992, 1816003).
Acknowledgments
Portions of this work were presented at the OSA Imaging and Applied Optics Congress in 2021, Compressive Spectral Imaging using Smoothness on Graphs [47]. Juan F. Florez thanks Fulbright Colombia and Colciencias for his doctoral fellowship, the University of Delaware’s graduate college for his dissertation fellowship, and Hoover Rueda, Carlos Mendoza, and Wenyi Ren for insightful and helpful technical discussions. This material is based upon work supported by the National Science Foundation under Grants NSF 1815992 and NSF 1816003.
Disclaimer. The images of the LEGO toys presented in the experimental results were captured by our imaging systems and were not endorsed by the LEGO trademark owners. Such images were used here as fair use to evaluate the spectral reconstruction capabilities of our method. LEGO is a trademark of the LEGO Group of companies, which does not sponsor, authorize or endorse the images in this paper. All Rights Reserved. https://www.lego.com/en-us/legal/notices-and-policies/fair-play/.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
Data availability
Data underlying the results presented in this paper are available in Ref. [48].
References
1. G. R. Arce, D. J. Brady, L. Carin, H. Arguello, and D. S. Kittle, “Compressive coded aperture spectral imaging: An introduction,” IEEE Signal Process. Mag. 31(1), 105–115 (2014). [CrossRef]
2. X. Yuan, D. J. Brady, and A. K. Katsaggelos, “Snapshot compressive imaging: Theory, algorithms, and applications,” IEEE Signal Process. Mag. 38(2), 65–88 (2021). [CrossRef]
3. J. Yang, X. Liao, X. Yuan, P. Llull, D. J. Brady, G. Sapiro, and L. Carin, “Compressive sensing by learning a Gaussian mixture model from measurements,” IEEE Trans. on Image Process. 24(1), 106–119 (2015). [CrossRef]
4. P. Meza, I. Ortiz, E. Vera, and J. Martinez, “Compressive hyperspectral imaging recovery by spatial-spectral non-local means regularization,” Opt. Express 26(6), 7043–7055 (2018). [CrossRef]
5. Y. Liu, X. Yuan, J. Suo, D. J. Brady, and Q. Dai, “Rank minimization for snapshot compressive imaging,” IEEE Trans. Pattern Anal. Mach. Intell. 41(12), 2990–3006 (2019). [CrossRef]
6. S. Zheng, Y. Liu, Z. Meng, M. Qiao, Z. Tong, X. Yang, S. Han, and X. Yuan, “Deep plug-and-play priors for spectral snapshot compressive imaging,” Photonics Res. 9(2), B18–B29 (2021). [CrossRef]
7. S. Arridge, P. Maass, O. Öktem, and C.-B. Schönlieb, “Solving inverse problems using data-driven models,” Acta Numerica 28, 1–174 (2019). [CrossRef]
8. V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,” IEEE Signal Process. Mag. 38(2), 18–44 (2021). [CrossRef]
9. Z. Meng, J. Ma, and X. Yuan, “End-to-end low cost compressive spectral imaging with spatial-spectral self-attention,” in European Conference on Computer Vision, pp. 187–204 (Springer, 2020).
10. D. Gilton, G. Ongie, and R. Willett, “Neumann Networks for Linear Inverse Problems in Imaging,” IEEE Trans. on Computational Imaging (2019).
11. G. Ongie, A. Jalal, C. A. Metzler, R. G. Baraniuk, A. G. Dimakis, and R. Willett, “Deep learning techniques for inverse problems in imaging,” IEEE J. Sel. Areas Inf. Theory 1(1), 39–56 (2020). [CrossRef]
12. L. Wang, T. Zhang, Y. Fu, and H. Huang, “Hyperreconnet: Joint coded aperture optimization and image reconstruction for compressive hyperspectral imaging,” IEEE Trans. on Image Process. 28(5), 2257–2270 (2019). [CrossRef]
13. 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, pp. 8032–8041 (2019).
14. D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” IEEE Signal Process. Mag. 30(3), 83–98 (2013). [CrossRef]
15. G. Puy, N. Tremblay, R. Gribonval, and P. Vandergheynst, “Random sampling of bandlimited signals on graphs,” Appl. Comp. Harmonic Anal. 44(2), 446–475 (2018). [CrossRef]
16. A. M. Tillmann, “On the computational intractability of exact and approximate dictionary learning,” IEEE Signal Process. Lett. 22(1), 45–49 (2015). [CrossRef]
17. B. Barazandeh and M. Razaviyayn, “On the behavior of the expectation-maximization algorithm for mixture models,” in 2018 IEEE Global Conference on Signal and Information Processing (GlobalSIP), pp. 61–65 (IEEE, 2018).
18. T. Glasmachers, “Limits of end-to-end learning,” in Asian Conference on Machine Learning, pp. 17–32 (PMLR, 2017).
19. M. E. Gehm, R. John, D. J. Brady, R. M. Willett, and T. J. Schulz, “Single-shot compressive spectral imaging with a dual-disperser architecture,” Opt. Express 15(21), 14013–14027 (2007). [CrossRef]
20. J. Monsalve, M. Marquez, I. Esnaola, and H. Arguello, “Compressive Covariance Matrix Estimation from a Dual-Dispersive Coded Aperture Spectral Imager,” in 2021 IEEE International Conference on Image Processing (ICIP), pp. 2823–2827 (IEEE, 2021).
21. X. Lin, Y. Liu, J. Wu, and Q. Dai, “Spatial-spectral encoded compressive hyperspectral imaging,” ACM Trans. Graph. 33(6), 1–11 (2014). [CrossRef]
22. E. Salazar, A. Parada-Mayorga, and G. R. Arce, “Spectral zooming and resolution limits of spatial spectral compressive spectral imagers,” IEEE Trans. Comput. Imaging 5(2), 165–179 (2019). [CrossRef]
23. A. Wagadarikar, R. John, R. Willett, and D. Brady, “Single disperser design for coded aperture snapshot spectral imaging,” Appl. Opt. 47(10), B44–B51 (2008). [CrossRef]
24. H. Arguello, C. V. Correa, and G. R. Arce, “Fast lapped block reconstructions in compressive spectral imaging,” Appl. Opt. 52(10), D32–D45 (2013). [CrossRef]
25. J. Zhang, D. Zhao, and W. Gao, “Group-based sparse representation for image restoration,” IEEE Trans. on Image Process. 23(8), 3336–3351 (2014). [CrossRef]
26. A. M. Teodoro, J. M. Bioucas-Dias, and M. A. Figueiredo, “A convergent image fusion algorithm using scene-adapted Gaussian-mixture-based denoising,” IEEE Trans. on Image Process. 28(1), 451–463 (2019). [CrossRef]
27. X. Dong, D. Thanou, M. Rabbat, and P. Frossard, “Learning graphs from data: A signal representation perspective,” IEEE Signal Process. Mag. 36(3), 44–63 (2019). [CrossRef]
28. G. Mateos, S. Segarra, A. G. Marques, and A. Ribeiro, “Connecting the dots: Identifying network structure via graph signal processing,” IEEE Signal Process. Mag. 36(3), 16–43 (2019). [CrossRef]
29. V. Kalofolias and N. Perraudin, “Large scale graph learning from smooth signals,” arXiv preprint arXiv:1710.05654 (2017).
30. V. Kalofolias, “How to learn a graph from smooth signals,” in Artificial Intelligence and Statistics, pp. 920–929 (2016).
31. X. Bresson and T. F. Chan, “Fast dual minimization of the vectorial total variation norm and applications to color image processing,” Inverse Probl. imaging 2(4), 455–484 (2008). [CrossRef]
32. M. F. Duarte and R. G. Baraniuk, “Kronecker compressive sensing,” IEEE Trans. on Image Process. 21(2), 494–504 (2012). [CrossRef]
33. C. V. Correa, H. Arguello, and G. R. Arce, “Snapshot colored compressive spectral imager,” J. Opt. Soc. Am. A 32(10), 1754–1763 (2015). [CrossRef]
34. H. Rueda-Chacon, F. Rojas, and H. Arguello, “Compressive spectral image fusion via a single aperture high throughput imaging system,” Sci. Rep. 11(1), 10311 (2021). [CrossRef]
35. A. Chakrabarti and T. Zickler, “Statistics of real-world hyperspectral images,” in CVPR 2011, pp. 193–200 (IEEE, 2011).
36. B. Arad and O. Ben-Shahar, “Sparse Recovery of Hyperspectral Signal from Natural RGB Images,” in European Conference on Computer Vision, pp. 19–34 (Springer, 2016).
37. 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 Trans. on Image Process. 20(3), 681–695 (2011). [CrossRef]
38. C. C. Paige and M. A. Saunders, “Solution of sparse indefinite systems of linear equations,” SIAM J. Numer. Anal. 12(4), 617–629 (1975). [CrossRef]
39. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. on Image Process. 13(4), 600–612 (2004). [CrossRef]
40. N. Yokoya, C. Grohnfeldt, and J. Chanussot, “Hyperspectral and multispectral data fusion: A comparative review of the recent literature,” IEEE Geosci. Remote Sens. Mag. 5(2), 29–56 (2017). [CrossRef]
41. L. Wang, Z. Xiong, D. Gao, G. Shi, and F. Wu, “Dual-camera design for coded aperture snapshot spectral imaging,” Appl. Opt. 54(4), 848–858 (2015). [CrossRef]
42. A. A. Wagadarikar, N. P. Pitsianis, X. Sun, and D. J. Brady, “Spectral image estimation for coded aperture snapshot spectral imagers,” in Image Reconstruction from Incomplete Data V, vol. 7076, p. 707602 (International Society for Optics and Photonics, 2008).
43. L. Wang, Z. Xiong, G. Shi, F. Wu, and W. Zeng, “Adaptive nonlocal sparse representation for dual-camera compressive hyperspectral imaging,” IEEE Trans. Pattern Anal. Mach. Intell. 39(10), 2104–2111 (2017). [CrossRef]
44. Y. Sun, Y. Yang, Q. Liu, and M. Kankanhalli, “Unsupervised Spatial-spectral Network Learning for Hyperspectral Compressive Snapshot Reconstruction,” IEEE Transactions on Geoscience and Remote Sensing (2021).
45. H. Dym, Linear algebra in action, vol. 78 (American Mathematical Soc., 2013).
46. X. Dong, D. Thanou, L. Toni, M. Bronstein, and P. Frossard, “Graph signal processing for machine learning: A review and new perspectives,” IEEE Signal Process. Mag. 37(6), 117–127 (2020). [CrossRef]
47. J. F. Florez-Ospina, D. L. Lau, K. Barner, and G. R. Arce, “Compressive Spectral Imaging using Smoothness on Graphs,” in Computational Optical Sensing and Imaging, pp. CTh2F–1 (Optical Society of America, 2021).
48. J. F. Florez-Ospina, A. K. M. Alrushud, D. L. Lau, and G. R. Arce, “Compressive spectral imaging using smoothness on graphs,” (2021). URL https://github.com/jfflorez/Compressive-spectral-imaging-using-graph-smoothness.git.
