Single-pixel imaging (SPI) is an emerging imaging technique employing a single-pixel detector, devoid of spatial resolution^[1]. SPI modulates the image of a scene using spatially resolved patterns to obtain a set of total light intensities, the correlation of which with the patterns yields the reconstructed images. The appeal of SPI lies in the efficiency of single-pixel detectors and post-processing algorithms, leading to its successful application in fields such as terahertz imaging^[2], fluorescence imaging^[3], and 3D imaging^[4]. Spectral imaging, another application of SPI, involves using a spectrometer as the single-pixel detector. Different from the conventional spectral imaging methods, spectral SPI acquires a 3D spectral cube using a set of patterns to modulate the image of an object, and then records the data by a single-pixel detector^[5]. This approach enables imaging with both high spectral and high spatial resolution. Despite its advantages, traditional spectral SPI requires long sequential measurements to acquire all the data, leading to long sampling time and extensive reconstruction period, particularly for high-spatial-resolution imaging.

The imaging efficiency of SPI can be enhanced through the incorporation of Hadamard basis reordering strategies and compressed sensing theory within sub-Nyquist SPI. By reordering a Hadamard pattern sequence, these strategies aim to reduce the required number of patterns, thereby shortening the SPI sampling time. However, the time consumption of applying compressed sensing algorithms for multiband image reconstruction at medium to high spatial resolution is excessively long, failing to meet the requirements for real-time imaging. To our knowledge, up to now, spectral SPI with spatial resolutions exceeding $128\times 128$ has not yet been achieved. Optimization strategies such as “Cakecutting,”^[6] “TV,”^[7] and “Zigzag,”^[8] which are simulated and experimented within the literature, typically restrict the spatial resolution to $128\times 128$ or lower. The overall sampling effectiveness of these strategies for higher resolutions and for multiband images remains unknown.

In response to these challenges, we propose an optimized Hadamard basis strategy. Initially, we arrange the Hadamard patterns in a 2D distribution based on a Gray code sequence. To avoid the cumbersome process inherent in traditional optimization strategies, we employ an index matrix to identify each pattern. Subsequently, we adopt a Morton scanning strategy to traverse the index matrix, directly generating the Hadamard pattern sequence. The unique characteristic of the Morton scanning curve allows spectral SPI to sample at a lower rate, significantly reducing the sampling time. For the task of multiband spectral reconstruction, we modify the traditional single-band fast Walsh-Hadamard transform (FWHT) to suit the multiband reconstruction tasks of spectral SPI and the forward measurement model based on Morton scanning. Ultimately, this results in an efficient reconstruction of high-resolution, multiband images in spectral SPI. The effectiveness of our proposed strategy is validated through simulations and experimental results, showing significant reduction in sampling and reconstruction time and demonstrating a promising solution for rapid and efficient spectral imaging.

As illustrated in Fig. 1, the spectral SPI system illuminates the scene with a sequence of modulating patterns while concurrently recording the total light intensity with a spectrometer. The spectral SPI measurement model is mathematically represented as $$\mathbf{Y}=\mathbf{A}\mathbf{X},$$(1)where $\mathbf{Y}\in {\mathbb{R}}^{M\times L}$ embodies the raw measurements obtained from the spectrometer. Here, $M$ denotes the number of patterns projected on DMD, while $L$ denotes the quantity of spectral channels. In this equation, $\mathbf{X}\in {\mathbb{R}}^{N\times L}$ represents the imaging spectral cube in a 2D form, $\mathbf{A}\in {\mathbb{R}}^{M\times N}$ signifies the measurement matrix, and $N$ indicates the total number of pixels of the image. The sampling ratio (SR) is thus defined as $\mathrm{SR}=M/N$. The conventional full sampling strategy utilizes the same number of patterns $M$ as the total number $N$ of pixels in an image. However, projecting all the required patterns can become prohibitively time-consuming if a high spatial resolution is desired. Therefore, a sub-sampling strategy ($M<N$) becomes a viable option for high-resolution spectral SPI.

Initially, we focus on the sampling process, which has the most significant impact on imaging efficiency. We utilize the widely used Hadamard matrix as our sampling matrix. Instead of constructing a Hadamard matrix, we analyze Hadamard basis vectors based on a Gray code sequence^[9]. First, let $V^{(k)}$ be the set of 1D Hadamard basis vectors of $k$-th order. $V^{(k)}$ can be constructed in a recursive way based on $V^{(k-1)}$. Specifically, each vector ${\mathbf{v}}^{(k-1)}\in V^{(k-1)}$ generates two basis vectors in $V^{(k-1)}$, which can be expressed as $$\begin{gathered} {\mathbf{v}}_m^{(k)}=[{\mathbf{v}}_n^{(k-1)},\alpha {\mathbf{v}}_n^{(k-1)}],n=1,\dots ,2^{k-1},m=1,\dots ,2^k, \\ \alpha \in \{-\mathrm{1,1}\}, \\ {\mathbf{v}}_n^{(k-1)}\in V^{(k-1)},{\mathbf{v}}_n^{(k)}\in V^{(k)}, \end{gathered}$$(2)where $[{\mathbf{v}}_1,{\mathbf{v}}_2]$ denotes the concatenation of ${\mathbf{v}}_1$ and ${\mathbf{v}}_2$. The initialization begins at $V^{(0)}=\{1\}$. The iterated construction of $V^{(k)}$ in Eq. (2) can also be illustrated using a binary tree. Figure 2(a) shows an example of the construction of set $V^{(2)}$. Each node at level $s$ represents a basis vector in $V^{(s)}$. All the nodes at depth $s$ compose the set $V^{(s)}$. Since $\alpha$ can be 1 or $-1$, each ${\mathbf{v}}_m^{(k-1)}$ in $V^{(k-1)}$ generates two unique ${\mathbf{v}}_n^{(k)}$ in $V^{(k)}$. It is observed from Fig. 2 that the branches of the tree are labeled with the values ${\alpha }_i$ to create the child basis vector so that the value of the child node can be generated by concatenating the father node and the multiplication of ${\alpha }_i$ and father node, as described in Eq. (2). Therefore, $[{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_k]$, called $\alpha$-index of ${\mathbf{v}}_i^k$, uniquely determines a basis vector in $V^{(k)}$.

It is noted that $V^{(k)}$ is an ordered set arranged according to the increasing order of the Gray code sequence, as shown in Fig. 2(b). Typically, Gray codes are represented by sequences of zeros and ones. Instead of using an $\alpha$-index to identify each pattern, a sequence $b=b_1{b}_2\dots b_k$ is employed. To streamline the construction of the sequence $b$, a direct method is used to determine the $b$ value of $v$ at the $i$-th position in the set $V^{(k)}$. This is achieved using the formula$$b=(i-1)\otimes ((i-1)\gg 1),$$(3)where $\otimes$ denotes the binary XOR operation, and $\gg$ signifies the logical right shift operation. For instance, as shown in Fig. 2(a), to calculate the $b$ corresponding to $[-1,-1]$ at the third position in $V^{(2)}$, one would apply the formula to compute the $10\otimes 01$, resulting in $b=11$.

We extend our consideration to two-dimensional Hadamard basis vectors. For the sake of brevity, all superscripts on variables are omitted in our discussion, with the understanding that the set and vectors considered are all under the $k$-th order. Let $P$ denote the 2D Hadamard basis set of $k$-th order. For any two elements ${\mathbf{v}}_i$ and ${\mathbf{v}}_j$ in $V^{(k)}$, their outer product ${\mathbf{v}}_i\otimes {\mathbf{v}}_j$ is an element of $P$. $P$ is composed of all possible outer products of pairs of elements from the set. The construction process of $P$ can be expressed as $$P=\{{\mathbf{v}}_i\otimes {\mathbf{v}}_j|{\mathbf{v}}_i,{\mathbf{v}}_j\in V^{(k)}\}.$$(4)

Figure 2(b) illuminates how to construct the set $P$ when $k=2$. For example, ${\mathbf{p}}_{12}=\left[\begin{array}{cccc}1& 1& -1& -1\\ 1& 1& -1& -1\\ 1& 1& -1& -1\\ 1& 1& -1& -1\end{array}\right]$ can be obtained by ${\mathbf{v}}_1\otimes {\mathbf{v}}_2$ when ${\mathbf{v}}_1=[\mathrm{1,1,1,1}]$ and ${\mathbf{v}}_2=[\mathrm{1,1,}\text{−}1,\text{−}1]$. Similar to the definition of $V^{(k)}$, the set of $P$ can be ordered in a 2D form as a $2^k\times 2^k$ matrix $\mathbf{p}$. Each element ${\mathbf{p}}_{ij}$ corresponds to a Hadamard pattern. Employing ${\mathbf{p}}_{ij}$ to represent the 2D Hadamard basis can be quite cumbersome. Here, we concatenate the Gray code sequence of $i$ and $j$, ${\mathbf{b}}_{ij}=[{\mathbf{b}}_i,{\mathbf{b}}_j],i,j=\mathrm{1,2,}\dots ,k$ to identify the ${\mathbf{p}}_{ij}$ at frequency coordinate $(i,j)$. The $2k$-length binary ${\mathbf{b}}_{ij}$ can be interpreted as a binary number of length $k$, which can be further converted into a decimal number, representing the location index $l$ of the pattern within a 1D natural order Hadamard pattern sequence^[10]. Therefore we construct an index matrix $\mathbf{I}=\left[\begin{array}{ccc}{I}_{11}& \dots & I_{1j}\\ ⋮& & ⋮\\ I_{i1}& \dots & I_{ij}\end{array}\right]$ to record location index $l$ of each ${\mathbf{p}}_{ij}$ of $P$. $I_{ij}$ can be computed as $$\begin{gathered} I_{ij}=\sum _{s=1}^{2k}{b}_s·2^{s-1}+1, \\ {\mathbf{b}}_{ij}=[b_1,b_2,b_s,\dots ,b_{2k}],s=\mathrm{1,2,}\dots ,2k. \end{gathered}$$(5)

For example, the ${\mathbf{b}}_{21}$ of the pattern ${\mathbf{p}}_{21}$ in position (2,1) is $[\mathrm{0,1,0,0}]$ and 3 is the row number of ${\mathbf{p}}_{21}$ in the natural order Hadamard matrix. The equation ${3=2}^0·0+2^1·1+2^2·0+2^3·0+1$ exists. When $k=2$, $\mathbf{I}=\left[\begin{array}{cccc}1& 3& 4& 2\\ 9& 11& 12& 10\\ 12& 15& 16& 14\\ 5& 7& 8& 6\end{array}\right]$. In Fig. 2(b), ${\mathbf{b}}_{ij}$ and $I_{ij}$ are shown below each ${\mathbf{p}}_{ij}$. This identification is crucial as we rearrange the detector sequence according to the natural order, enabling the use of a fast transform to accelerate the reconstruction of spectral data.

To achieve sub-sampling in spectral SPI, we aim to find an optimized pattern sequence within the 2D Hadamard pattern set $P$. We utilize the Morton scanning path^[11] to traverse $\mathbf{I}$ to obtain the set of location indices ${\mathrm{\Omega }}_M$ for $M$ optimized patterns. In Fig. 3, the Morton scanning path to traverse the $4\times 4$ Hadamard frequency domain is illustrated. Instead of traversing the 2D distributed Hadamard basis pattern set $P^{(2)}$, we traverse $\mathbf{I}$ to construct an ordered set ${\mathrm{\Omega }}_M$. For example, ${\mathrm{\Omega }}_4=\{\mathrm{1,3,9,11}\}$ when $M=4$. It is crucial to efficiently access and manipulate only the most significant frequency components. Morton scanning, by preserving spatial locality, allows for more efficient access to these critical components. Since the important frequency components (like low-frequency components in many images) are often spatially localized, Morton order ensures that these components are stored contiguously in SPI.

For SPI reconstruction under multi-band conditions, we employ a modified fast Walsh-Hadamard transform, termed HFWHT, to enhance the efficiency of spectral cube reconstruction. The details of HFWHT are shown in Fig. 4 and Algorithm 1. We initiate our process by applying zero-padding to the sub-sampling spectral data matrix ${\mathbf{Y}}_M$, adjusting its size from $M\times L$ to $N\times L$. ${\mathbf{Y}}_M$ is transformed to a measurement matrix ${\mathbf{Y}}_{nt}$ in natural order according to ${\mathrm{\Omega }}_M$. Subsequently, we perform $N/2$ addition and subtraction operations between the vectors in ${\mathbf{Y}}_{nt}$ similar to the process of FWHT. This step effectively processes each row’s vectors, containing information across various spectral bands simultaneously. Such an approach facilitates parallel computing in practical programming, significantly reducing the need for writing loops for spectral channels.

To evaluate the performance of the proposed strategy, we conducted numerical simulations using the CAVE spectral dataset^[13]. The dataset consists of 32 spectral images of indoor scenes, with a spectral cube resolution of $256\times 256\times 31$. The spectral channel ranges from 400 nm to 710 nm. We assessed the reconstruction quality of the spectral cube using two metrics: mean root mean squared error (MRMSE) and mean structural similarity index measure (MSSIM). These metrics were calculated as the mean values of the root mean squared error (RMSE) and the structural similarity index measure (SSIM)^[14] of each spectral channel.

First, we compared the proposed strategy with the “Zigzag” scanning path and two unoptimized orderings, “Natural” and “Random,” as well as two optimized orderings, “Cakecutting”^[7] and “Data”^[15]. While the original work predominantly performed the analysis at a resolution of 128, our simulations were conducted at a resolution of 256. To ensure the generalizability and broader applicability, simulations were carried out on 32 distinct cubes from the dataset. MRMSE and MSSIM are calculated for each reconstructed cube. We conducted such comparative simulations at intervals of 2% in the range from 1% to 50% sampling rates. Graphical representations of both MRMSE and MSSIM metrics are plotted to illustrate these comparisons in Fig. 5. Additionally, the standard deviation across multiple cubes was also depicted in the form of error bars in Fig. 5.

As observed in Fig. 5, Morton strategy consistently achieved the most optimal reconstruction results across various sampling rates in the simulations. “Zigzag” and “Cakecutting” emerged as relatively superior sampling strategies. Specifically, at a 25% sampling rate, “Zigzag” attained comparatively better reconstruction performance, while “Cakecutting” showed its relative strength at the same sampling rate. The performance of the “Data” was inferior, which may be attributed to the lack of strong similarity between the CAVE dataset and trained images. The unoptimized “Natural” and “Random” exhibited markedly poorer outcomes. Figure 6 presents a selection of the simulated reconstruction images from the aforementioned simulations, along with the corresponding MSSIM and MRMSE values. The figure displays pseudo-color images for three distinct cubes, reconstructed using four different sampling strategies at sampling rates of 2%, 10%, and 20%. The comparisons yield consistent conclusions. The Morton strategy obtains the best results, while “Zigzag” and “Data” achieve comparable outcomes, and “Natural” produces the least favorable results.

Next, we conduct the simulation of spectral cube reconstruction for a “balloons” cube under resolutions of $32\times 32$, $64\times 64$, $128\times 128$, and $256\times 256$. The sampling strategy consistently utilized is the Morton strategy. The computation time and MSSIM for different reconstruction strategies are shown in Tables 1 and 2. Reconstruction algorithms compared in this study encompass differential ghost imaging (DGI)^[16], Fourier domain regularization inversion (FDRI)^[17], and TVAL3^[18]. It is observed that the reconstruction times for DGI, FDRI, and TV surge significantly as the spatial resolution increases. The simulations are conducted under full sampling. It is necessary to reconstruct each of the 32 spectral bands separately and then sum up the reconstruction times. On the other hand, HFWHT, as a parallel reconstruction method, directly computes the complete cube. From Tables 1 and 2, we can conclude that HFWHT significantly outperforms DGI, FDRI, and TVAL3 in reconstruction time across various spatial resolutions. Meanwhile, the reconstruction performance of HFWHT is not inferior to traditional compressed sensing algorithms.

In this section, we validate the proposed strategy through experiment. The experimental setup, as illustrated in Fig. 7, includes a supercontinuum laser source (LEUKOS ROCK 400, 400–2400 nm), a collimating lens, a collection lens, a spectrometer (Ocean Optics FLAME-S-UV-VIS, 350–1000 nm), and a DMD (ViALUX, 9501). After collimation, the laser light illuminates the Rubik cube, and the light field is imaged through a double Gaussian structure TV lens of diameter 16 mm and focal length 50 mm, and is then modulated by the DMD. The total light is focused by the collection lens and recorded by the spectrometer. The spectral resolution of the spectral imaging system, contingent on the spectrometer’s slit width, is 2 nm in the range of 420–780 nm. The integration time of the spectrometer is set to 10 ms, considering the synchronization requirements of the spectrometer and DMD. Prior to the experiment, we utilized the proposed strategy to generate 4096 patterns with a resolution of $256\times 256$, reordered by the Morton scanning strategy, requiring a sampling time of 40.9 s. We obtained 181 spectral bands in the range of 420–780 nm. After processing the original spectrometer data, we obtained the measurement matrix ${\mathbf{Y}}_M$. Subsequently, we employ HFWHT and TVAL3 to reconstruct the target under a sampling rate of 5%. The retrieved spectral cube of 181 channels, spanning the range of 420 nm to 780 nm, is amalgamated into pseudo-color images as shown in Fig. 8(b) utilizing a color-matching method^[19]. The TVAL3 algorithm took 21 min, whereas HFWHT required only 0.41 s.

From the results displayed in Fig. 8, it is clearly evident that even at a 1% sampling rate, most of the original object’s image details can be reconstructed across all spectral bands. The low sampling rate undoubtedly leads to a significant increase in imaging efficiency. The primary reason for this is the suitability of our proposed Morton strategy for spectral SPI imaging scenarios at low sampling rates, as was prominently demonstrated in the simulations of the previous section. Next, we compare the effects of spectral reconstruction. The comparisons presented in Fig. 8 indicate that while the TVAL3 method achieves enhanced detail and reduced noise, HFWHT produces results that are not significantly inferior to TVAL3. This performance can be primarily attributed to the use of the Morton scanning strategy, whose properties are well-matched with the HFWHT reconstruction method. This advantage in reconstruction time becomes even more apparent with the increase in spectral bands and spatial resolution.

In conclusion, we propose a pipeline for efficient spectral SPI. To reduce sampling time, a strategy based on Morton scanning is incorporated within a 2D index matrix to produce a set of optimized Hadamard patterns at high spatial resolution. To facilitate efficiency in handling a substantial data volume, we leverage HFWHT to achieve spectral cube reconstruction, circumventing the need for handling extensive matrices. We demonstrate a spectral single-pixel imaging measurement system capable of capturing a $256\times 256\times 181$ spectral cube. When contrasted with traditional spectral single-pixel imaging methodologies, our imaging system demonstrates significantly lower sampling and reconstruction times. The sampling period is 40.9 s, while the reconstruction process is accomplished in 0.1 s. This reconstruction time represents an improvement of four orders of magnitude over conventional reconstruction methods. The sampling speed is predominantly constrained by the collection system, but can be increased through optimization of the data recording unit of the spectrometer. Furthermore, our strategy is applicable to other SPI modalities, 3D SPI for depth mapping, time-resolved SPI for fluorescence lifetime imaging, and so forth.