Spectral Zooming and Resolution Limits of Spatial Spectral Compressive Spectral Imagers

TL;DR

This paper develops a detailed model for the SSCSI, analyzing how aperture positioning affects spatial and spectral resolution, and demonstrates that physical aperture displacement enhances spectral detail in reconstructed datacubes.

Contribution

It introduces a rigorous discretization model for SSCSI, linking aperture displacement to spectral resolution improvements, supported by simulations and experiments.

Findings

A shift of the coded aperture increases spectral resolution.

The model accurately predicts resolution limits based on system parameters.

Experimental results validate the theoretical analysis.

Abstract

The recently introduced Spatial Spectral Compressive Spectral Imager (SSCSI) has been proposed as an alternative to carry out spatial and spectral coding using a binary on-off coded aperture. In SSCSI, the pixel pitch size of the coded aperture, as well as its location with respect to the detector array, play a critical role in the quality of image reconstruction. In this paper, a rigorous discretization model for this architecture is developed, based on a light propagation analysis across the imager. The attainable spatial and spectral resolution, and the various parameters affecting them, is derived through this process. Much like the displacement of zoom lens components leads to higher spatial resolution of a scene, a shift of the coded aperture in the SSCSI in reference to the detector leads to higher spectral resolution. This allows the recovery of spectrally detailed datacubes by physically displacing the mask towards the spectral plane. To prove the underlying concepts, computer simulations and experimental data are presented in this paper.

Figures

Tables

Table 1. TABLE I: Coherence of the sensing matrix for different values of s 𝑠 s

s	PSNR(dB)	$\mu \left(𝐀\right)$
0	26.7	1
0.01	28.53	0.997
0.02	30.22	0.971
0.03	30.35	0.971
0.05	30.71	0.922
0.07	31.54	0.888

Table 2. TABLE II: Theoretical and experimental spectral resolution for different s 𝑠 s parameters

s	${\Delta }_{\lambda }$, Theoretical (nm)	${\Delta }_{\lambda }$, Experimental (nm)
0.004	43	40
0.0078	22	28
0.011	16	20

Equations

$$R_1=\frac{bd}{b+d}\theta ,$$

$$\text{textsl}g(x,y)={\int }_{\Lambda }T(x(1-s)+s\alpha \lambda ,y)f(x,y,\lambda )d\lambda ,$$

$${\mathbf{g}}_{m,n}=\underset{x}{\int }\underset{y}{\int }{\int }_{\Lambda }T(x(1-s)+s\alpha \lambda ,y)f(x,y,\lambda )\times$$

$$\mathrm{rect}\left(\frac{x}{{\Delta }_d}-m,\frac{y}{{\Delta }_d}-n\right)d\lambda dydx,$$

$$\mathrm{rect}(x)=\left\{\begin{array}[]{ccc}1&\text{If}&0\leq x\leq 1\\ 0&&\text{otherwise,}\\ \end{array}\right.$$

$$T(x,y)=\sum _{m^{\prime }}\sum _{n^{\prime }}{\mathbf{t}}_{m^{\prime },n^{\prime }}\mathrm{rect}\left(\frac{x}{{\Delta }_c}-m^{\prime },\frac{y}{{\Delta }_c}-n^{\prime }\right),$$

$$\begin{gathered} {\mathbf{g}}_{m,n}=\underset{x}{\int }\underset{y}{\int }{\int }_{\Lambda }\sum _{m^{\prime }}\sum _{n^{\prime }}{\mathbf{t}}_{m^{\prime },n^{\prime }} \\ \times f(x,y,\lambda )\mathrm{rect}\left(\frac{x}{{\Delta }_d}-m,\frac{y}{{\Delta }_d}-n\right) \\ \times \mathrm{rect}\left(\frac{x(1-s)+s\alpha \lambda }{{\Delta }_c}-m^{\prime },\frac{y}{{\Delta }_c}-n^{\prime }\right)d\lambda dydx. \end{gathered}$$

$$\begin{gathered} {\mathbf{g}}_{m,n}=\underset{x}{\int }\underset{y}{\int }\underset{\Lambda }{\int }\sum _{m^{\prime }}\sum _{n^{\prime }}{\mathbf{t}}_{m^{\prime },n^{\prime }} \\ \times f(x,y,{\lambda }_o)\mathrm{rect}\left(\frac{x}{{\Delta }_d}-m,\frac{y}{{\Delta }_d}-n\right) \\ \times \mathrm{rect}\left(\frac{x(1-s)+s\alpha \lambda }{{\Delta }_c}-m^{\prime },\frac{y}{{\Delta }_c}-n^{\prime }\right)\delta (\lambda -{\lambda }_o)d\lambda dydx. \end{gathered}$$

$$\begin{gathered} {\mathbf{g}}_{m,n}=\underset{x}{\int }\underset{y}{\int }\sum _{m^{\prime }}\sum _{n^{\prime }}{\mathbf{t}}_{m^{\prime },n^{\prime }} \\ \times f(x,y,{\lambda }_o)\mathrm{rect}\left(\frac{x}{{\Delta }_d}-m,\frac{y}{{\Delta }_d}-n\right) \\ \times \mathrm{rect}\left(\frac{x(1-s)+s\alpha {\lambda }_o}{{\Delta }_c}-m^{\prime },\frac{y}{{\Delta }_c}-n^{\prime }\right)dydx. \end{gathered}$$

$$n\frac{{\Delta }_d}{{\Delta }_c}\le n^{\prime }\le (n+1)\frac{{\Delta }_d}{{\Delta }_c}-1,$$

$$\begin{gathered} \lfloor \frac{(m){\Delta }_d(1-s)+s\alpha {\lambda }_o}{{\Delta }_c}\rfloor \le m^{\prime }\le \\ \lfloor \frac{(m+1){\Delta }_d(1-s)+s\alpha {\lambda }_o}{{\Delta }_c}\rfloor , \end{gathered}$$

$${\mathbf{g}}_{m,n}=\sum _{m^{\prime }=m_l^{\prime }}^{m_r^{\prime }}\sum _{n^{\prime }=n\frac{{\Delta }_d}{{\Delta }_c}}^{(n+1)\frac{{\Delta }_d}{{\Delta }_c}-1}{W}_{m,m^{\prime }}\times {\mathbf{t}}_{m^{\prime },n^{\prime }}\times {\mathbf{f}}_{m^{\prime },n^{\prime },{\lambda }_o},$$

$$W_{m,m^{{}^{\prime}}}=\left\{\begin{array}[]{ccc}\frac{\frac{(m^{\prime}+1)\Delta_{c}-s\alpha\lambda_{o}}{1-s}-(m)\Delta_{d}}{\Delta_{c}/(1-s)}&\text{If}&m^{{}^{\prime}}=m_{l}^{{}^{\prime}}\\ \frac{(m+1)\Delta_{d}-\frac{m^{\prime}\Delta_{c}-s\alpha\lambda_{o}}{1-s}}{\Delta_{c}/(1-s)}&\text{If}&m^{{}^{\prime}}=m_{r}^{{}^{\prime}}\\ 1&\text{If}&m_{l}^{{}^{\prime}}<m^{{}^{\prime}}<m_{r}^{{}^{\prime}}\\ 0&&\text{otherwise}.\end{array}\right.$$

$$\begin{gathered} {\mathbf{g}}_{m,n}=\sum _{n^{\prime }=n\frac{{\Delta }_d}{{\Delta }_c}}^{(n+1)\frac{{\Delta }_d}{{\Delta }_c}-1}\left({\mathbf{t}}_{m^{\prime }-1,n^{\prime }}\times p_m+{\mathbf{t}}_{m^{\prime },n^{\prime }}\times (1-p_m)\right) \\ \times {\mathbf{f}}_{m,n^{\prime },{\lambda }_o}, \end{gathered}$$

$$p_{m}=\left\{\begin{array}[]{ccc}1&\text{If}&(m+1)\Delta_{d}\leq\frac{m^{\prime}\Delta_{c}-s\alpha\lambda_{o}}{(1-s)}\\ \frac{\frac{m^{\prime}\Delta_{c}-s\alpha\lambda_{o}}{(1-s)}-(m)\Delta_{d}}{\Delta_{d}}&\text{If}&(m+1)\Delta_{d}>\frac{m^{\prime}\Delta_{c}-s\alpha\lambda_{o}}{(1-s)}\end{array}\right.$$

$$\begin{gathered} {\mathbf{g}}_{m,n}=\underset{x}{\int }\underset{y}{\int }{\int }_{\Lambda }\sum _{m^{\prime }}\sum _{n^{\prime }}{\mathbf{t}}_{m^{\prime },n^{\prime }} \\ \times f(x_o,y_o,\lambda )\mathrm{rect}\left(\frac{x}{{\Delta }_d}-m,\frac{y}{{\Delta }_d}-n\right) \\ \times \mathrm{rect}\left(\frac{x(1-s)+s\alpha \lambda }{{\Delta }_c}-m^{\prime },\frac{y}{{\Delta }_c}-n^{\prime }\right) \\ \times \delta (x-x_o,y-y_o)d\lambda dydx, \end{gathered}$$

$$\begin{gathered} {\mathbf{g}}_{m,n}={\int }_{\Lambda }\sum _{m^{\prime }}\sum _{n^{\prime }}{\mathbf{t}}_{m^{\prime },n^{\prime }} \\ \times f(x_o,y_o,\lambda )\mathrm{rect}\left(\frac{x_o}{{\Delta }_d}-m,\frac{y_o}{{\Delta }_d}-n\right) \\ \times \mathrm{rect}\left(\frac{x_o(1-s)+s\alpha \lambda }{{\Delta }_c}-m^{\prime },\frac{y_o}{{\Delta }_c}-n^{\prime }\right)d\lambda , \end{gathered}$$

$$L=\lceil \frac{s\alpha ({\lambda }_{max}-{\lambda }_{min})}{{\Delta }_c}\rceil ,$$

$${\mathbf{g}}_{m,n}=\sum _{m^{\prime }=m_l^{\prime }}^{m_r^{\prime }}\sum _{n^{\prime }=n\frac{{\Delta }_d}{{\Delta }_c}}^{(n+1)\frac{{\Delta }_d}{{\Delta }_c}-1}\sum _{k=0}^{L-1}{W}_{m,m^{\prime }}{\mathbf{t}}_{m^{\prime }+k,n^{\prime }}{\mathbf{f}}_{m^{\prime },n^{\prime },k},$$

$$\begin{gathered} {\mathbf{g}}_{m,n}=\sum _{n^{\prime }=n\frac{{\Delta }_d}{{\Delta }_c}}^{(n+1)\frac{{\Delta }_d}{{\Delta }_c}-1}\sum _{k=0}^{L-1}({\mathbf{t}}_{m^{\prime }+k-1,n^{\prime }}\times p_m+ \\ {\mathbf{t}}_{m^{\prime }+k,n^{\prime }}\times (1-p_m)){\mathbf{f}}_{m,n^{\prime },k}, \end{gathered}$$

$$L=\lceil s\frac{N_d}{{\Delta }_c/{\Delta }_d}\beta \rceil ,$$

$${\vec{\mathbf{g}}}^{(q)}={\mathbf{H}}^{(q)}\mathbf{\Psi }\mathbf{\pi }q=1,2,...,Q,$$

$$\vec{\mathbf{g}}=\mathbf{H}\mathbf{\Psi }\mathbf{\pi },$$

$$\underset{\mathbf{\pi }}{\mathrm{argmin}}\left({\Vert \vec{\text{g}}-\text{A}\mathbf{\pi }\Vert }_2+\tau {\Vert \mathbf{\pi }\Vert }_1\right),$$

$$\begin{gathered} {\mathbf{g}}_{m,n}=\sum _{k=0}^{L-1}({\mathbf{t}}_{m^{\prime }+k-1,n^{\prime \prime }}\times p_m+ \\ {\mathbf{t}}_{m^{\prime }+k,n^{\prime \prime }}\times (1-p_m)){\mathbf{f}}_{m,n,k}, \end{gathered}$$

Full text

I Introduction

Hyperspectral imaging allows the simultaneous acquisition of spatial and spectral information of a given scene. It finds applications in a variety of fields such as remote sensing, medical imaging, food quality and spectroscopy [1, 2, 3]. As an alternative to the conventional hyperspectral sensing methods (push-broom, whisk-broom), the Coded Aperture Snapshot Spectral Imager (CASSI) [4, 5, 6, 7] was introduced; this is a Compressive Sensing based architecture [8], where the entire spatial-spectral datacube can be recovered from a few sensor measurements. A recently proposed enhancement of CASSI, replaces the binary on-off coded aperture by a pixelated color filter mask to achieve three dimensional coding prior to the optical sensing. This imager, known as colored CASSI [9, 10, 11, 12], increases both the quality of the reconstruction in the spatial domain and the attainable spectral resolution. Additionally, fewer snapshots are required to obtain accurate results [9]. The Snapshot Colored Compressive Spectral Imager (SCCSI)[13] and [14], on the other hand, represents a more compact version of the colored CASSI with remarkable performance, where the coded aperture is directly located on top of the sensor. Despite the advantages of both the colored CASSI and the SCCSI, the high fabrication cost of pixelated colored coded apertures limits their implementation in practice[15].

Recently, Lin et al. [16] proposed a new optical architecture called Spatial Spectral Compressive Spectral Imager (SSCSI). In this system, the spectral plane phenomena is exploited in order to carry out a three dimensional coding of the hyperspectral scene using a binary on-off mask. The spectral plane, see Fig. 1(a), is defined as the location where all the rays coming from the same wavelength meet each other in a single spatial point [17]. The SSCSI was first introduced as an attempt to replicate the functionality of the dual disperser CASSI, avoiding the calibration challenges, [18] and [19]. Coding in three dimensions, without the need of a colored mask, leads to high quality reconstructions while maintaining a low fabrication cost. Figure 1 depicts the ray tracing diagram of the SSCSI and the single disperser CASSI. The SSCSI system can be seen in Fig. 1(a), where the spatial-spectral coding process takes place after the spectral dispersion, and an in-focus image impinges on the sensor. In CASSI (see Fig. 1(b)), the coded aperture modulates an in-focus image of the scene and the dispersion occurs after the code modulation.

Some comparisons of the SSCSI with the different compressive spectral imagers have been done recently (Marín et al. in [20] and Medina-Rojas et al. in [21]). Nevertheless, a rigorous analysis of the SSCSI has not been developed, particularly on the resolution limits and the various parameters affecting the quality of the recovered scene. This paper provides such analysis with respect to the coded aperture and detector pitch sizes. Based on the continuous model studied in [16], a discrete expression of the sensor measurements is obtained. It is shown in section V-B and V-C, that spatial super-resolution is achievable, and a spectral zooming process can be done by displacing the coded aperture towards the spectral plane. Here, the term spectral zooming refers to the increment of the number of resolvable spectral bands on the recovered scene. This model is compared, by numerical simulations, against CASSI and colored CASSI in terms of spatial and spectral accuracy in the reconstruction. Additionally, a set of testbed experiments is performed in order to validate the underlying concepts. Some preliminary results of this research are described in [22] and [23].

This paper is organized as follows. In Section II, a description of the continuous model of SSCSI is presented. Section III presents the discretization of this model, offering a detailed description of the operators involved in the process. The structure of the sensing matrix is described in Section IV and the numerical simulations are presented in Section V. Sections VI and VII contain the experimental results and the conclusions, respectively.

II Sensing phenomena in SSCSI

This section contains the description of the sensing process for the SSCSI, developed in [16]. The SSCSI optical architecture depicted in Fig. 1(a), makes use of an objective lens to locate an image of the scene directly on the diffraction grating, where the spectral dispersion takes place. The rely lens is then used to get an in-focus image onto the sensor. Figure 2 depicts the ray propagation of the scene onto the spectral plane and into the sensor. Notice how all rays from a same wavelength converge to a single line on the spectral plane. The width of that spectral plane, $R_{1}$, is given by the following expression [17]

[TABLE]

where $b$ is the distance between the grating and the relay lens and $d$ is the distance between the spectral plane and the sensor (see Fig. 1(a)). The variable $\theta$ is the grating dispersion angle.

The parameter $s=d_{c}/d$ (see Fig. 1(a)), represents the position of the coded aperture with respect to the sensor. This normalized variable ($s=0$ when the coded aperture is on top of the sensor and $s=1$ when the coded aperture is on the spectral plane), strongly affects the coding process happening on the SSCSI. As depicted in Fig. 3(a), when $s=1$ each resolvable spectral band will be coded by a different column of the mask and a stretched version of the coded aperture will impinge on the sensor. When $0<s<1$ (see Fig. 3(b)) a similar process occurs, but in this case each resolvable band will be coded by a bigger portion of the coded aperture. When $s=0$, on the other hand (see Fig. 3(c)), all the spectral information will be coded by the same coded aperture pattern. As it will be described shortly, the parameter $s$ plays a crucial role in the resolution limits of SSCSI. The sensing process of the hyperspectral scene after being coded at a given $s$, is characterized by the following expression [16]

[TABLE]

where $\textsl{g}(x,y)$ is the continuous model representation of the information captured at the sensor (see Fig. 2), $T(x,y)$ is the coded aperture function, $f(x,y,\lambda)$ is the continuous representation of the hyperspectral scene, $\alpha$ is a parameter that indicates the position of a given wavelength on the spectral plane (see Appendix D for more details), and $\Lambda$ is the spectral range of interest, typically going from $400$nm to $700$nm [16]. Notice how the term $s$ appears on the argument of $T(x,y)$, which reflects the dependency of the coding process on that parameter. In this expression, it is assumed that the aperture of the objective lens is an infinitesimally small pinhole; however, in practice, an aperture with finite dimensions will introduce blurring effects as it will be explained in Section VI. Given that the sensor is composed of a discrete array of elements, a discretization process of Eq. (2) is a crucial step in understanding the physical phenomena in terms of the real captured information; this will be done in the following section.

III Discretization of the sensing phenomena

This section contains the steps followed in order to find a discrete expression for the captured data in the SSCSI. The discretization is performed based on the analysis of the product of the rectangular functions that represent both the coded aperture and the detector. The spatial resolution limits are first analyzed by assuming a single wavelength scene. This leads to a discrete measurement expression that is posteriorly extended to a multiple wavelength hyperspectral scene, by defining the spectral resolution limits on the SSCSI.

Consider the $(m,n)^{th}$ pixel detector, where $0\leq m,n\leq N_{d}-1$ are integer and unitless values, and $N_{d}\times N_{d}$ are the dimensions of the sensor array. The captured data on $(m,n)$ can be written as follows

[TABLE]

where $\Delta_{d}$ is the pitch size of a single detector element and $\mathrm{rect}(x)$ is a rectangular function defined as follows

[TABLE]

Likewise, the coded aperture can be written as

[TABLE]

where $\mathbf{t}_{m^{\prime},n^{\prime}}\in\{0,1\}\quad\forall\ (m^{\prime},n^{\prime})$, $0\leq m^{\prime},n^{\prime}\leq N_{c}-1$ are integer and unitless values, and $N_{c}\times N_{c}$ are the dimensions of the coded aperture mask. The variable $\Delta_{c}$ is the pitch size of the coded aperture. In this paper, the coded aperture full width is assumed to be the same as the detector full width, or $N_{c}\Delta_{c}=N_{d}\Delta_{d}$. Moreover, the detector and the mask are assumed to be totally aligned with respect to each other. Replacing (5) into (3) leads to

[TABLE]

A discretization measurement model of the last equation implies first a rigorous analysis of the spatial and spectral resolution limits related to the SSCSI, which will be done in the next subsections.

III-A Spatial Resolution Limits

Consider a spectral scene $f(x,y,\lambda_{o})\delta(\lambda-\lambda_{o})$ having just one spectral component, where $\delta(\lambda-\lambda_{o})$ is a Delta Dirac function located at $\lambda=\lambda_{o}$. Equation (6) can be rewritten as

[TABLE]

By solving the integral over $d\lambda$, Eq. (7) can be rewritten as

[TABLE]

The minimum spatially resolvable feature, given by Eq. (8) can be determined by establishing the intersecting region of the two rectangular functions in that expression. As it will be explained, this region is going to depend on the size of a coded aperture pixel projected onto the sensor.

III-A1 When $\frac{\Delta_{c}}{1-s}\leq\Delta_{d}$

Figure 4 left illustrates the coded aperture at a given position $s$ and its projection onto the sensor; the red lines delimit the detector elements. As depicted, a single coded aperture pixel of dimensions $\Delta_{c}\times\Delta_{c}$ is mapped into a $\frac{\Delta_{c}}{1-s}\times\Delta_{c}$ rectangle on the detector. This area defines the attainable spatial resolution of the recovered scene, as long as $\frac{\Delta_{c}}{1-s}\leq\Delta_{d}$. In this paper, it is assumed that $\Delta_{c}\leq\Delta_{d}$ and that the ratio $\Delta_{d}/\Delta_{c}$ is an integer greater than or equal to 1 (see Appendix E for $\Delta_{c}>\Delta_{d}$); this last condition implies that the mismatch between the coded aperture and the detector only applies to one axis (see Fig. 4, blue circled region).

Let $(m_{l}^{\prime},n^{\prime}_{l})$ and $(m^{\prime}_{r},n^{\prime}_{r})$ be the upper-leftmost and lower-rightmost coded aperture elements impinging on the $(m,n)^{th}$ detector element. The overlapping of the two rectangular functions in Eq. (8) can be seen in Figs. 5 and 6. From Fig. 5, it can be inferred that for the overlap to exist on the $y$ axis, the following condition must hold

[TABLE]

where the limits in the last inequality come from solving $n^{\prime}_{l}\Delta_{c}=n\Delta_{d}$ and $(n^{\prime}_{r}+1)\Delta_{c}=(n+1)\Delta_{d}$ for $n^{\prime}_{l}$ and $n^{\prime}_{r}$ respectively. Likewise, it can be shown based on Fig. 6, that for the two rectangular functions in Eq. (8) to overlap in the $x$ axis, the next condition must hold

[TABLE]

where $\lfloor\cdot\rfloor$ is the floor operator. The limits in the last inequality come from solving $\frac{m^{\prime}_{l}\Delta_{c}-s\alpha\lambda_{o}}{1-s}\leq m\Delta_{d}\leq\frac{(m^{\prime}_{l}+1)\Delta_{c}-s\alpha\lambda_{o}}{1-s}$ and $\frac{m^{\prime}_{r}\Delta_{c}-s\alpha\lambda_{o}}{1-s}\leq(m+1)\Delta_{d}\leq\frac{(m^{\prime}_{r}+1)\Delta_{c}-s\alpha\lambda_{o}}{1-s}$ for $m^{\prime}_{l}$ and $m^{\prime}_{r}$ respectively and considering that $m^{\prime}_{l}$ and $m^{\prime}_{r}$ are integer indexes. Given that the minimum resolvable spatial feature is of dimensions $\frac{\Delta_{c}}{(1-s)}\times\Delta_{c}$, and that the overlapping regions of the rectangular functions in Eq. (8) are defined in (9) and (10), the captured information in the $(m,n)^{th}$ detector element, $\mathbf{g}_{m,n}$, can be represented as

[TABLE]

where $m^{\prime}_{l}$ and $m^{\prime}_{r}$ are the limits of the expression in (10) and the variable $\mathbf{f}$ is the datacube representation of the hyperspectral scene, being $\mathbf{f}_{m^{\prime},n^{\prime},\lambda_{o}}$ the information of the datacube at position $(m^{\prime},n^{\prime})$ and spectral band $\lambda_{o}$; $\mathbf{f}_{m^{\prime},n^{\prime},\lambda_{o}}$ is obtained by applying the integral operators in Eq. (8). Notice that the spatial dimensions of $\mathbf{f}$, using an $N_{d}\times N_{d}$ sensor array, are defined as $\left\lceil N_{d}\frac{\Delta_{d}}{\Delta_{c}/(1-s)}\right\rceil\times N_{d}\frac{\Delta_{d}}{\Delta_{c}}$, where $\lceil\cdot\rceil$ is the ceiling operator. The parameter $W_{m,m^{\prime}}$ is the fraction of the coded aperture element $\mathbf{t}_{m^{\prime},n^{\prime}}$ impinging on the $(m,n)^{th}$ detector element, and it is defined as follows (see Appendix A for full derivation),

[TABLE]

A similar analysis was previously done for the CASSI system by Galvis et al. [24], in order to model the mismatching between the coded aperture and the sensor array. As an example, take the $(0,0)^{th}$ sensor element in Fig. 4 right; the captured information in this element will be equal to $\mathbf{g}_{0,0}=\sum_{m^{\prime}=0}^{2}\ \ \sum_{n^{\prime}=0}^{3}W_{0,m^{\prime}}\times\mathbf{t}_{m^{\prime},n^{\prime}}\times\mathbf{f}_{m^{\prime},n^{\prime},\lambda_{o}}$. In this particular case $W_{0,0}=1$ and $W_{0,1}=1$.

III-A2 When $\frac{\Delta_{c}}{1-s}>\Delta_{d}$

The projection of the coded aperture onto the sensor when $\frac{\Delta_{c}}{1-s}>\Delta_{d}$ is as shown in Fig. 7 left. Unlike the last scenario, the attainable spatial resolution in the $x$ axis is limited not by $\frac{\Delta_{c}}{1-s}$ but by the detector pitch size $\Delta_{d}$; the minimum spatially resolvable feature is then of dimensions $\Delta_{d}\times\Delta_{c}$. As when $\frac{\Delta_{c}}{1-s}\leq\Delta_{d}$ , a coded aperture-sensor mismatch is present. To deal with this particularity, a synthetic effective pattern is generated by the combination of elements proportional to the areas, as illustrated by the green circle in Fig. 7 right.

Given the $(m^{\prime},n^{\prime})^{th}$ coded aperture pixel and the $(m,n)^{th}$ detector element, the two rectangular functions in Eq. (8) intersect each other if $n^{\prime}$ falls on the interval specified in (9) and $m^{\prime}=\left\lfloor\frac{(m)\Delta_{d}(1-s)+s\alpha\lambda_{o}}{\Delta_{c}}\right\rfloor+1$ (see Appendix A). Taking into account that the minimum spatially resolvable feature is of dimensions $\Delta_{d}\times\Delta_{c}$, the captured information in a particular sensor element, $\mathbf{g}_{m,n}$, can be written as follows

[TABLE]

where $\mathbf{f}_{m,n^{\prime},\lambda_{o}}$ is obtained by applying the integral operators in Eq. (8). Notice that the spatial dimensions of $\mathbf{f}$, using an $N_{d}\times N_{d}$ sensor array, are defined as $N_{d}\times N_{d}\frac{\Delta_{d}}{\Delta_{c}}$. The expression $\tilde{\mathbf{t}}_{m,n^{\prime}}=\mathbf{t}_{m^{\prime}-1,n^{\prime}}\times p_{m}+\mathbf{t}_{m^{\prime},n^{\prime}}\times(1-p_{m})$ represents the effective coded aperture, where $p_{m}$ indicates the percentage occupied by $\mathbf{t}_{m^{\prime}-1,n^{\prime}}$ on the $(m,n)^{th}$ detector element and is defined as follows (see Appendix A),

[TABLE]

Unlike $W_{m,m^{\prime}}$, that represents a percentage value with respect to a coded aperture pixel, $p_{m}$ is a percentage value with respect to the detector pitch size $\Delta_{d}$. As an example, take the $(1,0)^{th}$ sensor element in Fig. 7 right. The captured information in this particular case will be equal to: $\mathbf{g}_{1,0}=\sum_{n^{\prime}=0}^{3}\left(\mathbf{t}_{0,n^{\prime}}\times p_{1}+\mathbf{t}_{1,n^{\prime}}\times(1-p_{1})\right)\times\mathbf{f}_{1,n^{\prime},\lambda_{o}}.$

Now that the spatial resolution limits were defined, an analysis to determine the minimum resolvable spectral band on the SSCSI must be developed in order to extend the discretization models, given in Eqs. (11) and (13), to the multi-spectral case.

III-B * Spectral Resolution Limits*

A similar process to the one done to calculate the spatial resolution must be done to find an expression for the spectral resolution in SSCSI. In this case, a multispectral scene $f(x_{o},y_{o},\lambda)\delta(x-x_{o},y-y_{o})$ representing a single spatial location, where $\delta(x-x_{o},y-y_{o})$ is a two dimensional Delta Dirac function located at $(x=x_{o},y=y_{o})$, is going to be considered. With this assumption, Equation (6) can be rewritten as follows

[TABLE]

By solving the integral over $dxdy$, Equation (15) can be rewritten as

[TABLE]

The attainable spectral resolution of the SSCSI is determined by the intersecting region of the two rectangular functions in Eq. (16). It can be shown that when $\Delta_{c}/(1-s)\leq\Delta_{d}$, that region has a length of $\Delta_{\lambda}=\frac{\Delta_{c}}{s\alpha}$ for certain values of $s$ (see Appendix B). Therefore, the maximum number of resolvable spectral bands is given by

[TABLE]

where $\lambda_{max}-\lambda_{min}=|\Lambda|$ is the spectral range of interest. Notice that the term $s\alpha(\lambda_{max}-\lambda_{min})$ indicates the spectral dispersion of the datacube for a given position $s$. Let $\lambda_{k+1}=\lambda_{min}+\frac{\Delta_{c}}{s\alpha}(k+1)$ and $\lambda_{k}=\lambda_{min}+\frac{\Delta_{c}}{s\alpha}(k)$ be two wavelengths, such that $\lambda_{k+1}-\lambda_{k}=\frac{\Delta_{c}}{s\alpha}$. The dispersion between these two bands is given by $s\alpha(\lambda_{k+1}-\lambda_{k})=\Delta_{c}$. This means that two adjacent and resolvable bands will be coded by the same coded aperture but shifted by one column $(\Delta_{c})$ as depicted in Fig. 8, top (for $k=0$). Eq. (11) can therefore be extended to the multispectral case as follows

[TABLE]

where $m^{\prime}_{l}=\left\lfloor\frac{(m)\Delta_{d}(1-s)+s\alpha\lambda_{min}}{\Delta_{c}}\right\rfloor$, $m^{\prime}_{r}=\left\lfloor\frac{(m+1)\Delta_{d}(1-s)+s\alpha\lambda_{min}}{\Delta_{c}}\right\rfloor$ and $W_{m,m^{\prime}}$ is defined in Eq. (12) (for $\lambda_{o}=\lambda_{min}$). The $k^{th}$ spectral band is defined here as the interval $[\lambda_{min}+\frac{\Delta_{c}}{s\alpha}k,\lambda_{min}+\frac{\Delta_{c}}{s\alpha}(k+1)]$; the shifting by one column between adjacent bands can be seen in the $``+k"$ term on the coded aperture array $\mathbf{t}$. The dimensions of the recovered datacube are given by $\left\lceil N_{d}\frac{\Delta_{d}}{\Delta_{c}/(1-s)}\right\rceil\times N_{d}\frac{\Delta_{d}}{\Delta_{c}}\times L$, where $N_{d}\times N_{d}$ are the dimensions of the sensor array and $L$ is given in Eq. (17). As an example, take the sensor element $(0,0)$ in Fig. 8, top; the captured information in that element related to the spectral band $[\lambda_{min},\lambda_{min}+\frac{\Delta_{c}}{s\alpha}]$ will be given by $\mathbf{g}_{0,0}=\sum_{m^{\prime}=0}^{2}\sum_{n^{\prime}=0}^{3}W_{0,m^{\prime}}\mathbf{t}_{m^{\prime},n^{\prime}}\mathbf{f}_{m^{\prime},n^{\prime},0}$. Likewise, the captured information for the spectral band $[\lambda_{min}+\frac{\Delta_{c}}{s\alpha},\lambda_{min}+\frac{2\Delta_{c}}{s\alpha}]$ will be given by $\mathbf{g}_{0,0}=\sum_{m^{\prime}=0}^{2}\sum_{n^{\prime}=0}^{3}W_{0,m^{\prime}}\mathbf{t}_{m^{\prime}+1,n^{\prime}}\mathbf{f}_{m^{\prime},n^{\prime},1}$. Notice that in this particular case, $W_{0,0}=1$ $W_{0,1}=1$. Since the actual shearing of the coded aperture seen at the sensor is approximated, this model can be called First order approximation. An analysis of the actual sheared coded aperture pattern and the parallel with the proposed model can be found in Appendix C.

When $\Delta_{c}/(1-s)>\Delta_{d}$, a similar analysis can be done to find that the attainable spectral resolution can be assumed as $\frac{\Delta_{c}}{s\alpha}$ (see Appendix B). Again, as when $\Delta_{c}/(1-s)\leq\Delta_{d}$, two adjacent and resolvable bands will be coded by the same coded aperture but shifted by one column, or $s\alpha(\lambda_{k+1}-\lambda_{k})=\Delta_{c}$, as depicted in Fig. 8, bottom (for $k=0$). Eq. (13) can therefore, be extended to the multispectral case as follows

[TABLE]

where the $k^{th}$ spectral band is defined as the interval $[\lambda_{min}+\frac{\Delta_{c}}{s\alpha}k,\lambda_{min}+\frac{\Delta_{c}}{s\alpha}(k+1)]$, $m^{\prime}=\left\lfloor\frac{(m)\Delta_{d}(1-s)+s\alpha\lambda_{min}}{\Delta_{c}}\right\rfloor+1$ and $p_{m}$ is as defined in Eq. (14) (for $\lambda_{o}=\lambda_{min}$). The shifting by one column between adjacent bands can be seen in the $``+k"$ term on the coded aperture array $\mathbf{t}$. The dimensions of the recovered datacube are given by $N_{d}\times N_{d}\frac{\Delta_{d}}{\Delta_{c}}\times L$, where $L$ is given in Eq. (17). As an example, take the sensor element $(1,0)$ in Fig. 8, bottom. The captured information in that element related to the spectral band $[\lambda_{min},\lambda_{min}+\frac{\Delta_{c}}{s\alpha}]$ will be given by $\mathbf{g}_{1,0}=\sum_{n^{\prime}=0}^{3}\left(\mathbf{t}_{0,n^{\prime}}\times p_{1}+\mathbf{t}_{1,n^{\prime}}\times(1-p_{1})\right)\mathbf{f}_{1,n^{\prime},0}$. Likewise, the captured information for the spectral band $[\lambda_{min}+\frac{\Delta_{c}}{s\alpha},\lambda_{min}+\frac{2\Delta_{c}}{s\alpha}]$ will be given by $\mathbf{g}_{1,0}=\sum_{n^{\prime}=0}^{3}\left(\mathbf{t}_{1,n^{\prime}}\times p_{1}+\mathbf{t}_{2,n^{\prime}}\times(1-p_{1})\right)\mathbf{f}_{1,n^{\prime},1}$. Again, the present model is an approximation and an analysis of the actual sheared coded aperture can be seen in Appendix C.

Figure 9 shows the coded aperture projected onto the sensor and the equivalent coded aperture for two adjacent bands when $\frac{\Delta_{c}}{1-s}>\Delta_{d}$ and $\frac{\Delta_{c}}{1-s}\leq\Delta_{d}$. Notice that the recovered datacube might not be square, as the resolution in each spatial axis can be different.

The dependency of the number of resolvable bands on $\Delta_{c}$ and $s$, as stated in Eq. (17), allows to enhance the spectral resolution of the system by tunning these parameters. Increasing $s$, which means moving the coded aperture towards the spectral plane, is equivalent to do a zooming process over the spectral dimension of the datacube. However, in experiments, placing the coded aperture at $s\neq 0$ decreases the spatial resolution of the recovered scene. This will be further discussed in the experimental part of the present paper.

An alternative expression for the number of the resolvable bands can be seen in the following equation (see Appendix D for full derivation),

[TABLE]

where $\beta=\alpha(\lambda_{max}-\lambda_{min})/(\Delta_{c}N_{c})$ is the ratio between the width of the spectral plane and the coded aperture full width. If $\beta=1$, the spectral plane fully occupies the coded aperture width. The variable $\beta$ is related to $R_{1}$, given by Eq. (1); therefore, a convenient way to increase $\beta$ is choosing a diffraction grating element with high diffraction angle $\theta$.

IV Matrix Forward Model

Given the discrete measurement model explained in the previous section, the sensing process can be expressed in a matrix form as follows

[TABLE]

where $Q$ is the number of captured shots and $\vec{\mathbf{g}}^{(q)}$ is the column-wise version of the measurements of length $N_{d}^{2}$. The variable $\boldsymbol{\pi}$ is a column-wise sparse representation of the datacube $\mathbf{f}$ in a basis $\boldsymbol{\Psi}$. If the datacube to be recovered has dimensions $N_{x}\times N_{y}\times L$, then the length of $\boldsymbol{\pi}$ is $N_{x}N_{y}L$. The structured matrix $\mathbf{H}^{(q)}$ performs the coding process of the datacube; its dimensions are $N_{d}^{2}\times N_{x}N_{y}L$. Figure 10 shows the structure of $\mathbf{H}^{(q)}$ when $\frac{\Delta_{c}}{1-s}\leq\Delta_{d}$ and $\frac{\Delta_{c}}{1-s}>\Delta_{d}$ for a single wavelength datacube. Each row of $\mathbf{H}^{(q)}$ corresponds to a sensor element, while each column corresponds to a datacube component. The values of the nonzero entries of the matrices are linked to Eqs. (18) and (19). Notice that one single sensor can measure several datacube elements.

For a multishot approach, the sensing process can be written as

[TABLE]

where $\mathbf{H}=\left[{\mathbf{H}^{(1)}}^{T},{\mathbf{H}^{(2)}}^{T},\dots,{\mathbf{H}^{(Q)}}^{T}\right]^{T}$ is a structured matrix with dimensions $QN_{d}^{2}\times N_{x}N_{y}L$, and $\vec{\mathbf{g}}=\left[\vec{\mathbf{g}}^{(1)T},\vec{\mathbf{g}}^{(2)T},\dots,\vec{\mathbf{g}}^{(Q)T}\right]^{T}$ is a vector of length $QN_{d}^{2}$. The structure of $\mathbf{H}$ for $Q=2$ and $L=2$ can be seen in Fig. 11. Unlike CASSI, the dispersion in the SSCSI is seen as the movement of the elements on the diagonals. An example on how to ensemble $\mathbf{H}$ when $\Delta_{d}=\Delta_{c}$ and the coded aperture has the same dimensions than the detector $(N_{d}=N_{c})$, can be seen in Algorithm 1. In this particular case, according to Eq. (9), $n^{\prime}=n$. The computational complexity of this algorithm is $\mathcal{O}(N^{2}LQ)$. When $\Delta_{c}/(1-s)<\Delta_{d}$, a similar process to the one specified in Algorithm 1 must be done. In this case, for every $(m,n)^{th}$ detector element, $n^{\prime}$, $m^{\prime}$ and $W_{m^{\prime},m}$ is calculated using Eqs. (9) and (10) and (12).

V Simulations

This section contains the simulation results in order to test the SSCSI first order approximation. The model is first compared to CASSI in terms of spatial and spectral quality in the reconstructions. After that, the spatial super-resolution and the spectral zooming concepts are tested in simulations. The performance of the SSCSI at different values of $s$ is then analyzed. A last subsection regarding the coherence of the sensing matrix and the impact on the quality of the reconstruction, is included. A hyperspectral scene was captured at the Computational Imaging and Spectroscopy Laboratory at University of Delaware, using a visible monochromator between 451nm and 642nm and a CCD camera, to posteriorly implement it as a ground-truth in simulations. The sparsity basis was chosen as $\boldsymbol{\Psi}=\boldsymbol{\Psi}_{DCT}\otimes\boldsymbol{\Psi}_{W}$, where $\boldsymbol{\Psi}_{W}$ is the 2D Wavelet Symlet 8 basis, $\boldsymbol{\Psi}_{DCT}$ is the Discrete Cosine basis and $\otimes$ is the Kronecker product [25]. The reconstruction is obtained as the solution of

[TABLE]

where $\left\|\cdot\right\|_{2}$ is the $\ell_{2}$ norm, $\left\|\cdot\right\|_{1}$ is the $\ell_{1}$ norm, $\tau$ is a regularization parameter that controls the sparsity of the solution and $\textbf{A}=\mathbf{H}\boldsymbol{\Psi}$. The Gradient Projection for Sparse Reconstruction algorithm (GPSR) is used in order to solve Eq. (23), as this algorithm has been shown to have optimal numerical performance and speed [26].

V-A Comparison to CASSI and Colored CASSI

The recovery process of a hyperspectral scene with $L=24$ bands was simulated with a $N_{d}\times N_{d}$ detector array, with $N_{d}=256$, and $\Delta_{d}=\Delta_{c}$ using the CASSI, the colored CASSI and the SSCSI first order approximation developed in this paper. Since $\Delta_{d}=\Delta_{c}$, $\Delta_{c}/(1-s)>\Delta_{d}$ for any value of $s$. Therefore, the spatial size of the recovered datacube is going to be $N_{x}\times N_{y}$, with $N_{x}=N_{y}=256$. Given that $L=24$ bands must be recovered, the coded aperture must be located at $s\approx 0.1$, according to Eq. (20) and assuming that $\beta=1$. The fact that $\beta=1$ means that $\lambda_{min}$ impinges in the leftmost side of the coded aperture and $\lambda_{max}$ impinges on the rightmost side of the coded aperture, or $\alpha(\lambda_{max}-\lambda_{min})=N_{c}\Delta_{c}$. The codes were generated such that the information captured in the shots is complementary, where the term complementarity is characterized by the condition $\sum_{q=1}^{Q}\mathbf{t}^{q}_{m,n}=1\ \forall\ m,n$; this structure, called Boolean, has been proven to be optimal for the CASSI architecture [27]; the transmittance of every coded aperture with boolean structure is inversely proportional to the number of captured shots. The regularization parameter $\tau$ was empirically chosen to attain optimal results for all the scenarios, being in the majority of cases between $\tau=1e-5$ and $\tau=1e-4$. The Peak Signal-to-Noise Ratio (PSNR) was the parameter implemented in order to determine the quality of the reconstructions. The PSNR is defined as $20\log_{10}\left(\max_{I}/\mathrm{MSE}^{1/2}\right)$, where $\max_{I}$ is the maximum possible value of the image and $\mathrm{MSE}$ is the mean squared error with respect to the ground-truth. Figure 12 shows the PSNR of the reconstructed scene as a function of the number of shots for the CASSI, the SSCSI first order approximation and the colored CASSI. The colored CASSI was implemented in two ways; in the first one, the coded aperture patterns were chosen such that the condition $\sum_{q=1}^{Q}\mathbf{t}^{q}_{m,n,k}=1\ \forall\ m,n,k$ holds, with no limitation on the type of optical filters; this is called ideal colored CASSI. In the second one, the types of filters are limited to four: low pass, high pass, band pass and band stop, and the patterns on the mask were chosen such that the condition $\sum_{q=1}^{Q}\mathbf{t}^{q}_{m,n,k}\geq 1\ \forall\ m,n,k$ holds. Here, $\mathbf{t}^{q}_{m,n,k}$ represents the $k^{th}$ component of the spectral response of the optical filter located at position $(m,n)$ and shot $q$. As depicted, the SSCSI exhibits much better performance than both, the CASSI and the colored CASSI with 4 filters while achieving similar results to the ideal colored CASSI. The comparison between the original scene and the recovered datacubes for $Q=3$ snapshots can be seen in Fig. 13, while a zoomed version can be observed in Fig. 14 left. Here, the hyperspectral scenes were mapped to the RGB domain for visualization purposes. A comparison between the original and recovered signatures for two different pixels, can be seen in Fig. 14 right. As depicted, the SSCSI, the ideal colored CASSI and the colored CASSI with 4 filters show spectrally accurate results.

V-B Spatial Super-Resolution

In [28], Arguello et al. proposed the reconstruction of spatially super-resolved hyperspectral scenes using CASSI, assuming that the ratio $\Delta_{d}/\Delta_{c}$ is an integer greater than 1. Likewise, with the SSCSI, Eq. (18) allows the reconstruction of spatially super-resolved datacubes with dimensions $\left\lceil N_{d}\frac{\Delta_{d}}{\Delta_{c}/(1-s)}\right\rceil\times N_{d}\frac{\Delta_{d}}{\Delta_{c}}\times L$. Consider $\Delta_{c}=\Delta_{d}/2$ and a $N_{d}\times N_{d}$ detector array, with $N_{d}=128$; with the SSCSI, a super-resolved datacube of spatial size $233\times 256$ and $L=24$ spectral bands, can be recovered if the coded aperture is located at $s\approx 0.1$ and $\beta=1$. Figure 15 top left and right, shows the recovered super-resolved datacube when the compression ratio is $CR=0.68$. The compression ratio, $CR$, is defined as the ratio between the number of captured measurements and the size of the reconstructed datacube, or $CR=\frac{QN_{d}^{2}}{N_{x}N_{y}L}$. Here, $N_{d}^{2}$ is the number of detector elements and $N_{x}\times N_{y}\times L$ is the size of the recovered datacube. The value of $CR$ must be less than $1$ in order to have compression during the capturing process. The $128\times 128\times 24$ non super-resolved datacube obtained from the same $128\times 128$ detector array can be seen in Fig. 15 bottom left and right. A simple inspection of the target reveals the effect of the super-resolution on the spatial details of the reconstruction. The PSNR of the reconstructed $233\times 256\times 24$ super-resolved datacube as a function of the compression ratio can be seen in Fig. 16 left. Here, the compression ratio is incremented by increasing the number of captured snapshots $Q$. A comparison between the original and the recovered spectral signatures can be seen in Fig. 16 right.

V-C Spectral Zooming

To illustrate the zooming process over the spectral dimension in a given scene, a spectral signature with strong peaks on three different wavelengths was introduced in the original datacube; the recovery was done for different values of $s$ and a compression ratio of $CR=0.2$. The conditions of $\beta=1$, $\Delta_{d}=\Delta_{c}$ and boolean coded apertures were assumed. As illustrated in Fig. 17, no spectral details can be identified when the recovery is done at $s\approx 0$ or $s\approx 0.1$; when $s=0.75$ and $s\approx 1$, on the other hand, the recovered spectral signatures are closer to the original one. This closeness is measured by the absolute value of the correlation coefficient $|r|$, between the original and reconstructed signatures.

V-D Performance at different $s$

As previously explained in section V-A, to recover a hyperspectral scene of $L=24$ bands using an $N_{d}\times N_{d}$ sensor, with $N_{d}=256$, and $\Delta_{d}=\Delta_{c}$, one must locate the coded aperture at $s\approx 0.1$, if $\beta=1$. This subsection analyses the impact of recovering that same scene for different values of $s$. Figure 18 shows the RGB profiles of the recovered scene for three different values of $s$ and $Q=6$ snapshots. It is clearly seen that with $s\approx 0$, the colors of the scene do not properly match the colors of the original datacube.

Figure 19 left, shows the PSNR of the reconstructed $256\times 256\times 24$ hyperspectral scene as a function of the coded aperture position, for $Q=6$ and $Q=12$ snapshots. As depicted, if $s\approx 0$, the quality of the reconstruction is low, since the system is not able to distinguish any spectral information of the scene. Increasing $s$ leads to enhance the quality of the reconstruction. However, when $s$ gets close to $1$, the PSNR drops as it was also shown in [16]. The comparison between the spectral signatures for different values of $s$ can be seen in Fig. 19 right; reducing $s$ to $0.02$ makes the fine details of the curve to be lost while maintaining its shape; for $s\approx 0$, no spectral information is recovered.

V-E Coherence of the sensing matrix

In compressive sensing, the coherence of the sensing matrix $\mu(\mathbf{A})$, where $\mathbf{A}=\mathbf{H}\boldsymbol{\Psi}$, has been extensively used as a measure of the quality of the sensing process, where matrices with a low coherence lead to high quality reconstructions [29, 30, 31, 32]. Table I contains the PSNR values of the reconstructions and the coherence of the sensing matrix $\mu(\mathbf{A})$ for different values of $s$. Here, a $N_{d}\times N_{d}$ sensor array with $N_{d}=64$ was considered, with $\Delta_{d}=\Delta_{c}$ and $\beta=1$. A multispectral scene of dimensions $64\times 64\times 8$ was recovered and the coded apertures exhibit boolean structure. As it can be appreciated in Table I, $\mu(\mathbf{A})$ is connected to the quality of the recovered scene and low values of $\mu(\mathbf{A})$ are associated to high PSNR values. Further analysis of the coherence in SSCSI will be done in future work.

VI Experimental measurements

The SSCSI was experimentally implemented as depicted in Fig. 20. A TAMROM AF 70-300mm lens locates an image of the scene at the transmissive diffraction grating (300 grooves/mm). The diffracted image is then reshaped using a 4f system composed of two 75mm, 2” lenses. In a 4f system, the image located at a distance $f$ with respect to the first lens is reshaped at a distance $f$ with respect to the second lens. The 4f system provides major flexibility in the movement of the coded aperture, in the sense that $s=0$ is accessible, since it is located in a spatial position far from the sensor. The spectral plane will be located in between the lenses that compose the 4f system. The coded scene is then reshaped at the sensor by a 1”, 35mm lens. The implemented sensor is a StingrayTM $640\times 480$ CCD monochrome camera with 9.9µm pitch size. The spectral range of interest was defined as 480nm-620nm.

VI-A Spectral Resolution

The theoretical spectral resolution can be determined by defining the size of the spectral plane with respect to the coded aperture ($\beta$ in (20)). It was found that the 480nm-620nm spectral range occupies a physical space three times the coded aperture width, therefore, $\beta\approx 3$. This was done by locating a white board on the spectral plane and measuring the physical distance between $\lambda_{min}=480$nm and $\lambda_{max}=620$nm. The number of resolvable bands with $\Delta_{c}/\Delta_{d}=1$ and $N_{d}=N_{c}=256$ is given by $L=\left\lceil 256s\beta\right\rceil$. On the other hand, the criteria used to determine whether two adjacent spectral bands are resolvable in experiments, can be seen in Fig. 21. As depicted, the patterns that code two adjacent spectral bands at a given $s$ are first captured with the sensor; then, the $i^{th}$ rows from those patterns are compared. In order for the two bands to be resolvable, the curves must be shifted by one pixel. Table II shows the experimental and theoretical spectral resolution for different values of $s$.

VI-B Experimental Datacube Reconstruction

A hyperspectral scene was reconstructed using $2$ snapshots and different values of $s$. The coded aperture patterns exhibit boolean structure and were printed on a photo-mask with pitch size of $\Delta_{c}=19.4$µm. The matrix $\mathbf{H}$ was generated by first locating a white board target and then capturing the coded aperture pattern for each spectral band to be recovered, at a given $s$. Then, that pattern is located in the respective diagonal of the matrix $\mathbf{H}$; this is done for every captured snapshot. After that, the white board is replaced by the hyperspectral scene, and the snapshots are properly captured. The reconstruction algorithm is posteriorly executed to find the optimal coefficients $\boldsymbol{\pi}$. The RGB profiles of the reference scene and the reconstructed datacubes can be seen in Fig. 22, while the spectral bands and recovered spectral signatures are contained in Fig. 23 and Fig. 24. As depicted, increasing $s$ allows the spectral information of the scene to be better distinguished.

The SSCSI bases its functionality on the location of the coded aperture in an out of focus position or $s>0$, as depicted in Fig. 1. As mentioned in Section II, the present model assumes an infinitesimally small aperture in the objective lens. However, in experiments, the finite size of the aperture introduces blurring effects when $s>0$ and, in fact, increasing $s$ deteriorates the spatial quality of the reconstructed hyperspectral scene. Figure 22 shows the zooming over one portion of the recovered scene for the different values of $s$. It is evident that the colors of the scene become more noticeable for a bigger $s$, but at the same time, the spatial quality drops.

There are several factors to consider in the experiments. Besides the blurring proportional to $s$, that affects the spatial quality of the reconstructions, one might find other challenges when implementing the SSCSI. One of them is the amount of light that passes through the system, since the aperture of the objective lens must be as small as possible. At the same time, the size of this aperture directly affects the resolution of the SSCSI, given that the minimum resolvable feature for any optical system is determined by the diffraction limit. On the other hand, the fact that the spectral plane is three times bigger than the coded aperture $(\beta\approx 3)$, leads to extreme wavelengths to impinge out of the mask; this means, in practice, that the range of wavelengths is reduced to the ones impinging inside the mask.

VII Conclusions and further steps

This papers develops a rigorous discretization measurement model of the Spatial Spectral Compressive Spectral Imager (SSCSI), based on the physical dimensions of the coded aperture and the detector array, and the spectral dispersion introduced by the diffraction grating, characterized by $\alpha$. Two different scenarios were proposed:

• •

$\Delta_{d}\leq\Delta_{c}/(1-s)$: This scenario is modeled by Eq. (18), where $W_{m,m^{\prime}}$ is used to model the mismatch between the coded aperture and the detector array. The spatial resolution is defined as $\Delta_{c}/(1-s)\times\Delta_{c}$. As described in Section V-B, a spatially super-resolved datacube of dimensions$\left\lceil N_{d}\frac{\Delta_{d}}{\Delta_{c}/(1-s)}\right\rceil\times N_{d}\frac{\Delta_{d}}{\Delta_{c}}\times L$, can be recovered when this condition holds.

• •

$\Delta_{d}>\Delta_{c}/(1-s)$: This scenario is modeled by Eq. (19), where an effective mask is generated by the combination of coded aperture pixels proportional to the areas, using the variable $p_{m}$. The spatial resolution is defined as $\Delta_{d}\times\Delta_{c}$

The number of spectral bands that can be recovered using the SSCSI was found to depend on three variables: The coded aperture pitch size $\Delta_{c}$, the spectral dispersion introduced by the diffraction grating characterized by $\alpha$, and the coded aperture position represented by the parameter $s$. In fact, it was found that a bigger $s$ increases the number of resolvable spectral bands. Therefore, a movement of the coded aperture towards the spectral plane can be seen as performing a spectral zooming over a given hyperspectral scene. The SSCSI performance, in experiments, depends on the aperture of the objective lens. Increasing the size of the aperture makes the spatial quality of the recovered scene to drop as the coded aperture is displaced towards the spectral plane. On the other hand, reducing the aperture also decreases the amount of light passing through the system.

Ongoing work includes the optimization of the coded aperture patterns for SSCSI based on the concept of coherence and the proposed model in this paper. Similarly, a more accurate discretization model is being analyzed.

Acknowledgment

This research project was funded by the Department of Homeland Security, Science and Technology Directorate (Contract $\mathrm{\#}$HSHQDC-15-C-B0016).

Appendix A

$W_{m,m^{\prime}}$ and $\tilde{\mathbf{t}}_{m,n^{\prime}}$ full calculation

$W_{m,m^{\prime}}$

As mentioned in Section III-A, $W_{m,m^{\prime}}$ indicates the fraction of the coded aperture $\mathbf{t}_{m^{\prime},n^{\prime}}$ impinging on the sensor element $(m,n)$. A graphical description of $W_{m,m^{\prime}}$ can be seen in Fig. 25.

As depicted, three cases can be distinguished; when $m^{\prime}=m^{\prime}_{l}$ (gray area in Fig. 25), $W_{m,m^{\prime}_{l}}$ can be calculated as $W_{m,m^{\prime}_{l}}=\frac{\frac{(m^{\prime}_{l}+1)\Delta_{c}-s\alpha\lambda_{o}}{1-s}-(m)\Delta_{d}}{\Delta_{c}/(1-s)}$. The term $\Delta_{c}/(1-s)$ is used as a normalization parameter that allows one to obtain the percentage of $\mathbf{t}_{m^{\prime}_{l},n^{\prime}}$ measured by the sensor element $(m,n)$. Likewise, when $m=m^{\prime}_{r}$ (yellow region in Fig. 25), $W_{m,m_{r}^{\prime}}$ can be calculated as $W_{m,m^{\prime}_{r}}=\frac{(m+1)\Delta_{d}-\frac{m^{\prime}_{r}\Delta_{c}-s\alpha\lambda_{o}}{1-s}}{\Delta_{c}/(1-s)}$. Finally, if $m^{\prime}_{l}<m^{\prime}<m^{\prime}_{r}$ (green area in Fig. 25), $W_{m,m^{\prime}}$ can be calculated as $W_{m,m^{\prime}}=\frac{\Delta_{c}/(1-s)}{\Delta_{c}/(1-s)}=1$. The fact $W_{m,m^{\prime}}=1$ indicates that $\mathbf{t}_{m^{\prime},n^{\prime}}$ is fully impinging on sensor element $(m,n)$. Equation (12) contains a general expression for $W_{m,m^{\prime}}$ considering all the mentioned cases. An example of how to calculate $W_{m,m^{\prime}}$ can be done based on Fig. 4 right. Here, $W_{0,2}=\frac{\Delta_{d}-\frac{2\Delta_{c}}{(1-s)}}{\Delta_{c}/(1-s)}$, $W_{1,2}=\frac{\frac{3\Delta_{c}}{(1-s)}-\Delta_{d}}{\Delta_{c}/(1-s)}$ and $W_{0,1}=1$.

$\tilde{\mathbf{t}}_{m,n^{\prime}}$

The effective coded aperture $\tilde{\mathbf{t}}_{m,n^{\prime}}$ is given by the expression $\tilde{\mathbf{t}}_{m,n^{\prime}}=\mathbf{t}_{m^{\prime}-1,n^{\prime}}\times p_{m}+\mathbf{t}_{m^{\prime},n^{\prime}}\times(1-p_{m})$, where $n^{\prime}$ is given by Eq. (9), $m^{\prime}=\left\lfloor\frac{(m)\Delta_{d}(1-s)+s\alpha\lambda_{o}}{\Delta_{c}}\right\rfloor+1$ and the parameter $p_{m}$ indicates the percentage of the $(m,n)^{th}$ sensor element occupied by the $(m^{\prime}-1,n^{\prime})^{th}$ coded aperture element. Figure 26 contains a graphical explanation to calculate $p_{m}$. Notice that, according to this figure, $m\Delta_{d}\leq\frac{m^{\prime}\Delta_{c}-s\alpha\lambda_{o}}{1-s}\leq(m+1)\Delta_{d}$, or equivalently $m^{\prime}=\left\lfloor\frac{(m)\Delta_{d}(1-s)+s\alpha\lambda_{o}}{\Delta_{c}}\right\rfloor+1$ (considering that $m^{\prime}$ is an integer index).

The portion of the $(m,n)^{th}$ sensor element occupied by $\mathbf{t}_{m^{\prime}-1,n^{\prime}}$ (gray area in Fig. 26) can be calculated as $p_{m}=\frac{\frac{m^{\prime}\Delta_{c}-s\alpha\lambda_{o}}{(1-s)}-(m)\Delta_{d}}{\Delta_{d}}$, where the sensor pitch size $\Delta_{d}$ appears in the denominator for normalization purposes. If $\frac{m^{\prime}\Delta_{c}-s\alpha\lambda_{o}}{1-s}\geq(m+1)\Delta_{d}$, $\mathbf{t}_{m^{\prime}-1,n^{\prime}}$ fully occupies the $(m,n)^{th}$ sensor element and therefore $p_{m}=1$. A general expression to calculate $p_{m}$, including the two cases mentioned can be seen in Eq. (14).

Appendix B

Spectral Resolution calculation

$\frac{\Delta_{c}}{1-s}\leq\Delta_{d}$

The spectral resolution of SSCSI is given by the length of the region where Eq. (16) is different from zero. This is defined by the overlapping of the two rectangular functions in that expression, which can be seen in Fig. 27. As depicted two concrete cases can be analyzed. In the first one (gray area in Fig. 27), $m\Delta_{d}\leq\frac{(m^{\prime}_{1}+1)\Delta_{c}-s\alpha\lambda}{1-s}$ and $m\Delta_{d}\geq\frac{m^{\prime}_{1}\Delta_{c}-s\alpha\lambda}{1-s}$ or, in other words, $\frac{m^{\prime}_{1}\Delta_{c}-m\Delta_{d}(1-s)}{s\alpha}\leq\lambda\leq\frac{(m^{\prime}_{1}+1)\Delta_{c}-m\Delta_{d}(1-s)}{s\alpha}$. This interval has an extension of $\Delta_{\lambda}=\frac{\Delta_{c}}{s\alpha}$. Notice that the same $\Delta_{\lambda}$ is obtained for the yellow region in Fig. 27.

In the second case (green area), $m\Delta_{d}\leq\frac{m^{\prime}_{2}\Delta_{c}-s\alpha\lambda}{1-s}$ and $(m+1)\Delta_{d}\geq\frac{(m^{\prime}_{2}+1)\Delta_{c}-s\alpha\lambda}{1-s}$, or equivalently $\frac{(m^{\prime}_{2}+1)\Delta_{c}-(m+1)\Delta_{d}(1-s)}{s\alpha}\leq\lambda\leq\frac{m^{\prime}_{2}\Delta_{c}-m\Delta_{d}(1-s)}{s\alpha}$, which has an extension of $\Delta_{\lambda}=\frac{\Delta_{d}(1-s)-\Delta_{c}}{s\alpha}$.

Let the ratio $\Delta_{d}/\Delta_{c}=C$, where $C$ is an integer greater than 1. In order to make $\frac{\Delta_{c}}{s\alpha}>\frac{\Delta_{d}(1-s)-\Delta_{c}}{s\alpha}$, one can prove that $s>1-\frac{2}{C}$. For $C=2$ (the case analyzed in section V-B), the inequality $\frac{\Delta_{c}}{s\alpha}>\frac{\Delta_{d}(1-s)-\Delta_{c}}{s\alpha}$ holds for $s>0$. Therefore, $\Delta_{\lambda}=\frac{\Delta_{c}}{s\alpha}$ can be taken as an upper bound of the spectral resolution. When $C>2$, $\frac{\Delta_{c}}{s\alpha}<\frac{\Delta_{d}(1-s)-\Delta_{c}}{s\alpha}$ for certain values of $s$, and $\Delta_{\lambda}=\frac{\Delta_{c}}{s\alpha}$ represents the exact spectral resolution given by the SSCSI.

$\frac{\Delta_{c}}{1-s}>\Delta_{d}$

A similar analysis to the previous one can be done to find the spectral resolution in this case. Figure 28 shows the overlapping of the two rectangular functions given in Eq. (16). Again two scenarios must be analyzed, In the first one (gray area in Fig. 28), $\frac{(m^{\prime}+1)\Delta_{c}-s\alpha\lambda}{1-s}\leq(m+1)\Delta_{d}$ and $\frac{(m^{\prime}+1)\Delta_{c}-s\alpha\lambda}{1-s}\geq(m)\Delta_{d}$, or $\frac{(m^{\prime}+1)\Delta_{c}-(m+1)\Delta_{d}(1-s)}{s\alpha}\leq\lambda\leq\\ \frac{(m^{\prime}+1)\Delta_{c}-(m)\Delta_{d}(1-s)}{s\alpha}$, which has an extension of $\Delta_{\lambda}=\frac{\Delta_{d}(1-s)}{s\alpha}$. In the second case (gray area in Fig. 28 occupies the whole sensor element), $\frac{(m^{\prime})\Delta_{c}-s\alpha\lambda}{1-s}\leq m\Delta_{d}$ and $\frac{(m^{\prime}+1)\Delta_{c}-s\alpha\lambda}{1-s}\geq(m+1)\Delta_{d}$, or equivalently $\frac{m^{\prime}\Delta_{c}-m\Delta_{d}(1-s)}{s\alpha}\leq\lambda\leq\frac{(m^{\prime}+1)\Delta_{c}-(m+1)\Delta_{d}(1-s)}{s\alpha}$, which has an extension of $\Delta_{\lambda}=\frac{\Delta_{c}-\Delta_{d}(1-s)}{s\alpha}$. Notice that, since $\frac{\Delta_{c}}{1-s}>\Delta_{d}$, $\frac{\Delta_{c}}{s\alpha}>\frac{\Delta_{d}(1-s)}{s\alpha}$ and $\frac{\Delta_{c}}{s\alpha}>\frac{\Delta_{c}-\Delta_{d}(1-s)}{s\alpha}$. Therefore, $\frac{\Delta_{c}}{s\alpha}$ can be taken as an upper bound of the spectral resolution.

Appendix C

Actual coded aperture shearing at the sensor

Let recall Eq. (6), which characterizes the SSCSI sensing process. By analyzing the $x$ dimension of the coded aperture rectangular function in this expression, given by $\mathrm{rect}\left(\frac{x(1-s)+s\alpha\lambda}{\Delta_{c}}-m^{\prime},\frac{y}{\Delta_{c}}-n^{\prime}\right)$, one can obtain a figure that depicts the actual sheared pattern. For a given $m^{\prime}$, $x$ is limited by $\frac{m^{\prime}\Delta_{c}-s\alpha\lambda}{1-s}\leq x\leq\frac{(m^{\prime}+1)\Delta_{c}-s\alpha\lambda}{1-s}$. Figure 29(a) shows the sheared coded aperture seen at the sensor for both cases, when $\Delta_{c}/(1-s)\leq\Delta_{d}$ and $\Delta_{c}/(1-s)>\Delta_{d}$, in the $x-\lambda$ plane. Notice how $m^{\prime}$ is bounded by the limits given above. The dashed red lines in this figure represent the spectral division every $\frac{\Delta_{c}}{s\alpha}$, while the dashed black lines indicate the limits of the $m^{th}$ sensor element. The first order approximation, given by Eqs. (18) and (19), can be seen in Fig. 29(b). Notice how, for example, the spectral band $[\lambda_{1},\lambda_{2}]$ is coded by the same pattern that codes the hyperspectral scene at $\lambda_{1}$.

Appendix D

Derivation of Eq. (20)

The number of resolvable bands in the SSCSI system is given by $L=\frac{s\alpha(\lambda_{max}-\lambda_{min})}{\Delta_{c}}$, where $\alpha(\lambda_{max}-\lambda_{min})$ is the width of the spectral plane. Let $\beta=\alpha(\lambda_{max}-\lambda_{min})/N_{c}\Delta_{c}$, where $N_{c}\Delta_{c}$ is the coded aperture width. The number of resolvable bands as a function of $\beta$ can be rewritten as $L=s\beta N_{c}$. If it is assumed that $N_{d}\Delta_{d}=N_{c}\Delta_{c}$, which means that the coded aperture width is the same as that of the sensor width, last expression can be written as Eq. (20).

The ceiling operator $\lceil\cdot\rceil$ in Eqs. (17) and (20) was introduced in order fully use the spectral resolution given by the SSCSI.

Appendix E

Discretization Model for $\Delta_{c}>\Delta_{d}$

This paper proposes a SSCSI discretization measurement model for the specific case of $\Delta_{d}/\Delta_{c}=C$, where $C$ is an integer greater than or equal to 1. If $\Delta_{c}=C_{2}\Delta_{d}$, where $C_{2}$ is an integer greater than 1, the condition $\Delta_{c}/(1-s)>\Delta_{d}$ holds for any value of $s$ and Eq. (19) can be rewritten as follows:

[TABLE]

where $m^{\prime}=\left\lfloor\frac{(m)\Delta_{d}(1-s)+s\alpha\lambda_{min}}{\Delta_{c}}\right\rfloor+1$, $p_{m}$ is defined in Eq. (14) (for $\lambda_{o}=\lambda_{min}$), and $n^{\prime\prime}=\left\lfloor\frac{n}{C_{2}}\right\rfloor$. Notice that the minimum spatially resolvable element in the reconstructed scene is given by $\Delta_{d}\times\Delta_{d}$, and the spatial size of the recovered datacube is equal to $N_{d}\times N_{d}$.