1.

Introduction

Hyperspectral imaging (HSI) refers to the acquisition of multiple images corresponding to more spectral bands than just the typical three color channels red, green, and blue in RGB imaging. Early applications of HSI included determining the spectra of stars and their composition¹ or the identification of clay minerals from airborne measurements.² During the last decade HSI has been becoming more important in areas such as medicine,³ food inspection,⁴ environmental monitoring,⁵ and many others. Multiple approaches for hyperspectral data acquisition have been proposed in the literature and the reader might refer to an extended historical overview of HSI in Ref. 6 and different HSI systems in Ref. 7.

Computed tomography imaging spectrometer (CTIS) has been first proposed by Okamoto and Yamaguchi⁸ and independently by Bulygin and Vishnyakov.⁹ A three-dimensional (3D) data cube ($x$, $y$, and wavelength) is determined on the basis of projections taken from different “directions.” The advantages of CTIS are its hyperspectral snapshot capability and simple setup.¹⁰ The disadvantages are a large unused sensor area and a time-consuming iterative algorithm that is required to reconstruct the hyperspectral image.

2.

Computed Tomography Imaging Spectrometer

A standard CTIS design is shown in Fig. 1. The object is illuminated by a broadband light source. A field stop is used to block the outer part of the intermediate image and thereby determines the systems’ field of view. For our work, we use an aperture with a quadratic shape which results in a quadratic hyperspectral image. Two further lenses first collimate the light coming from the field stop and then image it onto the sensor. In between is a dispersive element that separates the light depending on the wavelength. Commonly, a diffractive optical element (DOE) (e.g., made of several crossed linear gratings) is employed.¹² We use a binary computer-generated hologram (CGH) that leads to an image shown in Fig. 2(a). In the central part of the image, the original broadband scene is located. Outside of this central area, we have 12 dispersed copies due to the first diffraction order of the CGH. It can be seen that the longer wavelengths, i.e., red region, are diffracted more than the shorter ones, i.e., blue region. This way the spectral information is spatially encoded on the sensor. Since we use a monochromatic sensor we obtain an image as shown in Fig. 2(b). The CGH is calculated using the direct binary search algorithm.¹³

Fig. 1

Optical schematic of a standard CTIS setup. An intermediate image of an object located to the left is created and cropped by a field stop (the butterfly image shown here is taken from Ref. 11). The cropped image is then diffracted into several directions by the DOE and re-imaged onto the sensor on the right. The CTIS sensor image is shown in Fig. 2.

Fig. 2

CTIS sensor image: (a) RGB representation (butterfly image taken from Ref. 11) is generated using the methods described in Sec. 3.2.1. The colors of the zeroth-order are changed compared with the ground truth since the spectral efficiencies of the DOE are taken into account. (b) monochromatic sensor image as used to extract the hyperspectral information. Reconstructions from this image are shown in Sec. 4. The diffraction orders are slightly brightened here for illustration purposes.

Research on CTIS has focused on reducing the unused sensor space with different diffraction patterns,^14–17 accelerating the reconstruction process,^18–23 and finding applications for CTIS in, e.g., fluorescence microscopy,²⁴ infrared spectroscopy,^25,26 and astronomy.^27,28

3.

Neural Networks for CTIS

In CTIS, the spatial information (zeroth diffraction order) and the spectral information (higher diffraction orders) are projected to different locations on the image sensor whereas large areas remain unused. These circumstances should be considered when designing a neural network for processing CTIS images.

Douarre et al.³⁰ took into account the characteristics of a CTIS image and separated the zeroth diffraction order from the higher diffraction orders. The higher diffraction orders are cropped and rotated to align them next to each other. The realigned diffractions are then processed by rectangular convolution kernels. The zeroth diffraction order is processed independently from the higher ones and is concatenated later. However, by rotating the higher diffraction orders, the underlying spatial information is also rotated and the neural network had the burden to learn to connect the right areas from each diffraction order.

Huang et al.²⁹ cropped the different diffraction orders and stacked them on top of each other as a new dimension. Since they only worked with horizontal and vertical diffractions, five subimages were cropped. The cropping eliminates most of the empty sensor space and no spatial information is altered due to rotation. However, the corresponding spectral information for one pixel is still far apart on the different layers, which requires additional and unnecessary work for the neural network and could be avoided by applying prior system knowledge to the CTIS image.

3.1.

Proposed Architecture

In our proposed projection to cube (p2cube) architecture, the higher diffraction orders are reshaped into cubes before convolutions are applied. The reshaping is done in a way that classic convolution can connect areas where the corresponding spectral information is distributed. The idea is adapted from deformable convolutions where the kernels can connect different regions of the input.³¹ While the deformation in deformable convolutions is learned during backpropagation we apply a fixed deformation since the pattern is known from the PSF. Even though this reshaping is rather simple, it plays a pivotal role in a robust reconstruction of the hyperspectral cube from the CTIS image.

The p2cube architecture is shown in Fig. 3 with a brief overview of the applied reshaping in the red box. Here, the input image is reshaped into 13 cubes. One cube for every higher diffraction order and one cube, where the zeroth diffraction order is replicated to fit the spectral dimensions of the other cubes and it provides the spatial information. The reshaped cubes are then processed first by a 3D deconvolution layer and then four convolutional layers which together form the baseline architecture of our approach. We added a U-Net³² like extension shown in the gray box and the combined architecture will hereafter be referred to as p2cubeU. The number of trainable parameters from the two architectures is shown in Table 1. In comparison, the architecture from Huang et al.²⁹ has over 85 million parameters used to reconstruct hyperspectral images with a comparable size and the architecture from Douarre et al.³⁰ uses around 16 million parameters for a classification task on CTIS images.

Fig. 3

Schematic visualization of the p2cube architecture with the reshaping layer visualized in more detail in the red box, and the U-Net like extension in the gray box, the combination of both architectures is referred to as p2cubeU.

Table 1

Trainable parameters of the neural networks and number of iterations for the EM algorithm with the resulting computation time. The EM algorithm is running on an Nvidia RTX 2070 SUPER and the neural networks are tested on an Nvidia GTX 1060 6 GB.

	EM	p2cube	p2cubeU
Parameters	—	575,291	1,566,091
Iterations	200	—	—
Computation time	11 s	0.14 s	0.17 s

The reshaping is visualized in more detail in Fig. 4 for two higher diffraction orders. Several squares are cropped and stacked on top of each other. The step size between the squares can be controlled in the model. A step size of five pixels shows good results in combination with a five-by-five kernel size. On the right side of Fig. 4, the reshaped cubes are shown. The white squares visualize the area which is connected by a convolution kernel. The connected area from the cube is projected back to the original diffraction order. It can be seen that the classic convolution, which is applied to the cube, connects exactly the area where the spectral information from a pixel in the object space is distributed along. This way even the diagonal diffraction orders can be processed by rectangular convolution kernels. By applying a 3D deconvolution to all thirteen cubes, all areas where the spectral information is distributed are connected.

Fig. 4

p2cube reshaping in more detail: only five layers are shown for better visualization. The area connected by a convolution kernel corresponds to the area in which the spectral information is distributed. The size of the kernel is exaggerated for better visualization.

3.2.

Data Generation

The training data consists of sensor images as they are captured with a CTIS system together with the hyperspectral representation of the scenes (ground truth). In the following, we describe the synthetic data generation pipeline and our CTIS system that was used to acquire a real dataset with sensor images and corresponding hyperspectral ground truth.

3.2.1.

CTIS-Emulation for synthetic data generation

We created a synthetic data set by emulating the CTIS system using Fourier optics to calculate the PSF for every wavelength (color) channel as described in Ref. 33.

A monochromatic point source is considered as a scene in the far-field, which results in a uniform wavefront in the Fourier space. The CGH introduces a spatially dependent phase shift according to the optical path length. Our binary CGH creates a phase shift of $\pi$ at 550 nm. The image in the Fourier space is cropped differently for every wavelength to keep the pixel size in the image plane constant and, therefore, allows the simulation of the CTIS image without interpolation in the spatial dimension. We use a sensor resolution of $1200\times 1200\text{pixels}$. An inverse Fourier transformation leads to the wavefront in the image plane. By taking the square of the complex magnitude from the wavefront, the intensity at the sensor (PSF) is obtained. This has to be done for every wavelength. Afterward, each wavelength-dependent PSF is convolved with the corresponding spectral band from the hyperspectral cube. By summing up all spectral channels the final grayscale sensor image is obtained. We use a spectral resolution of one nanometer to create the PSF. The hyperspectral images are interpolated in the spectral dimension to fit the calculated PSF. It has to be noted that this method does not take optical aberrations into account and assumes a spatially invariant system. Only diffraction due to the finite lens diameter is considered.

Hyperspectral images from three different publicly available datasets are used: BGU iCVL Hyperspectral Image Dataset,³⁴ CAVE Multispectral Image Database,³⁵ and TokyoTech 31-band Hyperspectral Image Dataset.¹¹ We extracted the bands from 450 to 690 nm with a spectral resolution of 10 nm, which leads to 25 spectral channels. The original images are first randomly split into training (76.5%), validation (19%), and test (4.5%) sets. Then, using a sliding window, nine subimages were extracted from a given image and then resized to a spatial resolution of $100\times 100\text{pixels}$. A tenth image is the whole original one but resized to $100\times 100\text{pixels}$ also. This way the amount of samples used to train the network is increased by a factor of ten to a total number of 2585 hyperspectral images.

3.2.2.

Real world CTIS data

Even though the emulation of the CTIS system using Fourier optics takes into account effects such as wavelength dependent diffraction efficiency, it does not include all optical aberrations from the lenses and internal reflections in the setup. To test the performance of the proposed neural network for real-world data, a real dataset of CTIS images paired with hyperspectral ground truth cubes was captured. A schematic of the expanded optical setup is shown in Fig. 5. A second optical path is added by inserting a beam splitter after the field stop. There, the intermediate image is again collimated and spectrally filtered by a tunable color filter.³⁶ A second sensor records the spectral channels sequentially. All images were recorded from 455 to 695 nm with a step size of 10 nm resulting in 25 spectral channels. The CTIS images have a resolution of $1088\times 1088\text{pixels}$ and the hyperspectral ground truth images were downsampled to match the resolution of the zeroth diffraction order of $139\times 139\text{pixels}$. The recorded real-world dataset consists of 495 images that were split into 400 images for training, 80 for validation, and 15 images for testing.

Fig. 5

Schematic of our setup to capture CTIS images and corresponding hyperspectral cubes. By inserting a beam splitter after the field stop, both sensors capture the intermediate image during the field stop. The tunable color filter allows to record the spectral channels from 455 nm to 695 nm with a 10-nm step size.³⁶

4.

Results

In the following, only results that are part of the test data set and thus not utilized during training or validation are discussed. In addition, to assess the influence of noise, we use one data set without noise and another data set consisting of the same images with shot noise added to the detector images. The neural network is trained with and without noise accordingly. For the EM algorithm, a noise-free calibration (PSF image) is used for both sets. The number of iterations for the EM algorithm was fixed to 200 since the performance stagnated afterward. Hyperspectral reconstruction from real CTIS images with the EM algorithm was not successful and is therefore not presented. The failure of the EM algorithm is due to a higher noise level in the sensor images and the requirement of a spatially invariant PSF, which is not fulfilled by the real CTIS setup.

To evaluate the reconstruction performance, we used three different quality criteria as presented in Tables 2 and 3. An individually optimized scaling factor is applied to the reconstructed images from the EM algorithm to deliver the best result. This is necessary since the sensor images are normalized and the total energy is lost. This is valid because only relative and not absolute values are of interest anyway. As the first quantitative criterion, we used the root mean squared error (RMSE). Since the data range is between 0 and 1, the error can be directly interpreted. Similar to the RMSE is the mean absolute error (MAE), which, in contrast to the RMSE, takes individual strong deviations (peaks) less into account. The mean structural similarity index measure (MeanSSIM) on the other hand is a completely different metric that is usually used to determine the quality of color images. In our case, we use the mean values over all channels. This, therefore, does not account for spectral characteristics anymore.⁴⁰ However, we consider it is still a valuable metric since it is different than the others. Respective results are shown in Table 2.

Table 2

Quantitative comparison between EM, p2cube, and p2cubeU architectures for synthetic data. The values labeled with noise are derived from the noisy sensor images. For RMSE and MAE a small value close to 0 shows good agreement of the reconstruction with the ground truth data. For MeanSSIM, in contrast, a value close to 1 shows a good structural similarity of the data.

Table 3

Quantitative results on real CTIS images using p2cube and p2cubeU architectures. Reconstruction with EM failed in this case since the PSF is not spatially invariant and the CTIS image is too noisy.

It can be seen that p2cubeU delivers the best results for all metrics for noise-free and noisy images. The EM algorithm is only slightly worse for the noise-free case. For noisy images, however, the quality of the images reconstructed with the EM algorithm decreases significantly. Here, the p2cube and p2cubeU architectures show their strength and can deliver results that are close to the noise-free ground truth images.

Qualitative comparisons are shown in Figs. 6 and 7. The RGB representation of the different images is obtained via the CIE 1931 filter response functions. Notice that the reconstruction from noise-free CTIS images works pretty well with all three methods. However, all images show some imperfections. For the EM reconstruction, one can see straight artifacts that correspond to the direction of diffraction of the DOE. This is typical and can often be observed with such kind of algorithms. With the neural network, it is sometimes noticeable that the color of small details is not reproduced correctly, e.g., the small stripes in the lower left part of the butterfly wing. For the emulated images with noise, it can be seen that the neural network-based approach delivers better results since the images from the EM algorithm become very noisy in the spectral dimension.

Fig. 6

RGB representation of the ground truth and the reconstructed hyperspectral images from the synthetic noise-free sensor images. The images are $100\times 100\text{pixels}$ and have 25 channels. Spectra of two points per image are shown on the right.

Fig. 7

RGB representation of the ground truth and the reconstructed hyperspectral images from the synthetic noisy sensor images. The images are $100\times 100\text{pixels}$ and have 25 channels. Spectra of two points per image are shown on the right.

In Fig. 8, reconstructed samples from the real CTIS data are shown. Notice that the images from p2cube and p2cubeU show very few artifacts in the spatial dimension. Furthermore, p2cubeU achieves better performance and produces well-reconstructed images inline with the corresponding ground truth.

Fig. 8

RGB representation of the ground truth and the reconstructed hyperspectral images from the real-world images. The images are $139\times 139\text{pixels}$ and have 25 channels.

5.

Conclusion and Outlook

In this paper, we presented a deep convolutional neural network for hyperspectral image reconstruction from synthetic as well as real CTIS data. Since spectral and spatial information overlap on the sensor image, it is rearranged in a first step so that the subsequent convolutional layers can better link spatially related areas. This reshaping organizes the CTIS image in a more efficient way for the network than others have demonstrated so far and it ensures a robust training/reconstruction with over 30 times less parameters than Douarre et al.³⁰ for their classification task and 100 times less parameters than Huang et al.²⁹ for their hyperspectral image reconstruction. Our approach produces high-quality reconstruction results on synthetic and real CTIS data.

The advantages of the proposed neural network are the reduced computation time for inference by a factor $>60$ compared with the conventional EM algorithm and the reduced complexity compared with other networks dealing with CTIS data.^29,30 With a computation time between 0.14 and 0.17 s per image on a GTX 1060 6 GB, it enables real-time performance with common hardware.

In our experiments, the spatial resolution of the hyperspectral images is the same as the spatial resolution of the zeroth diffraction order. However, it should be possible to generate higher spatial resolution images from CTIS images since the spatial information is also multiplexed in the higher diffraction orders. This will be the focus of future work.