Open Access
Add to BibSonomy
Abstract
We present a fast reconstruction algorithm for hyperspectral images, utilizing a small amount of data without the need for any training. The method is implemented with a dual disperser hyperspectral imager and makes use of spatial-spectral correlations by a so-called separability assumption that assumes that the image is made of regions of homogenous spectra. The reconstruction algorithm is simple and ready-to-use and does not require any prior knowledge of the scene. A simple proof-of-principle experiment is performed, demonstrating that only a small number of acquisitions are required, and the resulting compressed data-cube is reconstructed near instantaneously.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
Corrections
Elizabeth Hemsley, Ibrahim Ardi, Tony Rouvier, Simon Lacroix, Hervé Carfantan, and Antoine Monmayrant, "Fast reconstruction of hyperspectral images from coded acquisitions using a separability assumption: erratum," Opt. Express 31, 43521-43521 (2023)https://opg.optica.org/oe/abstract.cfm?uri=oe-31-26-43521
1. Introduction
Hyperspectral (HS) imagers provide precise spectral information for every pixel in a scene, and are useful for a range of applications - e.g. deep space or earth observation, material analysis, gas detection, etc. The time taken to obtain the 3D HS data-cube can be long, particularly for scanning type imagers, necessitating several acquisitions with a 2D sensor, and the volume of data produced is large. Often much of the information in the HS datacube is redundant, given that a real world scene contains a considerable amount of spectral-spatial correlations [1,2]. These correlations lead to the field of compressed sensing, which reduces the amount of data needed to reconstruct the HS datacube, based on the assumption of an underlying sparse dictionary [3–7]. More recently, neural networks have been trained to reconstruct a HS cube using a similarly small amount of data [8–10]. But these methods either require significant time during the reconstruction, or a learning phase which make them not straightforward and ready-to-use. In this paper we present an alternative method requiring no training or prerequisite knowledge of the scene. The amount of data required for this algorithm is similarly small, and the reconstruction time is very short. A fundamental point is that we do not attempt to reconstruct the entire hyperspectral datacube, that is to reconstruct the spectrum for each pixel of the scene. The key idea of our approach is to observe that nearby pixels often share the same spectrum with only a small panchromatic intensity difference, for example due to slightly different lighting conditions. If two pixels have differing spectra, the difference in the panchromatic intensity will likely be more significant. Therefore the scene can be divided into spectrally homogenous regions using information from the panchromatic image, and we need only assign a spectrum to each region. This approach has two main advantages over more widespread compressed sensing approaches; firstly the reduction in computing cost, and secondly the possibility to independently reconstruct regions within the scene, leading to a reduced size in the resulting datacube. Some approaches using panchromatic or RGB side information have been proposed to improve accuracy for compressed sensing schemes [11,12], but this often requires increased complexity of optical setup and additional calibration challenges. Particularly for single disperser systems, an additional arm and camera must be incorporated and accurately co-registered with the rest of the data, which can present a significant challenge. The advantage of a system which is perfectly co-registered by design, such as the DD-CASSI, is that the panchromatic corresponds spatially to the rest of the data, measuring the sum over the wavelength bands, with no additional calibration. These advantages allow us to define this particular reconstruction algorithm, which is much more computationally efficient than classical compressed sensing algorithms.
There are other systems which have the potential to fulfill this quality of co-registration, e.g. the more recently developed SSCSI system [13,14] or the CASSI-SLM [11]. However, the SSCSI system requires a very small input aperture [15], so light collection efficiency is problematic. For the CASSI-SLM system, obtaining the panchromatic as a sum over the wavelength bands is similarly challenging.
Paper outline: the next section depicts how the reconstruction problem for a dual disperser HS imager is given by a least square minimization, augmented with a spectral regularization. Section 3 presents the experimental setup, how the scene ground truth is recovered, and the metric used to asses the performance of our approach. Section 4 analyzes experimental results for one simple and one more complex scene, and in Section 5 we discuss how to remove reliance on segmentation of the panchromatic image.
2. Methodology
Due to spectral-spatial correlations within a HS scene, adjacent pixels are likely to share the same spectra. Following this logic, we can separate any hyperspectral scene into spectrally homogenous regions, defined as spatially connected areas where the spectra for every pixel is the same, apart from a multiplication factor which depends on the intensity. The whole homogenous region can thus be described by one single spectrum and an intensity map given by the panchromatic intensity. We denote this assumption as the separability assumption (SA). If a homogenous region in the scene consists of $K$ pixels indexed by $k \in \{1 \ldots K\}$, the HS information $\mathbf {o}_k$ of pixel $k$ in the region is given by:
where $\mathbf {s}$ is the normalized spectrum of the region made of $W$ spectral bands, $P_k$ the panchromatic intensity of that pixel, and the index k corresponds to all the pixels of the region. Therefore, the hyperspectral information $\textbf {o}$ of the whole region can be written:The measurement is implemented using a dual disperser HS imager [16], illustrated in Fig. 1. For this architecture each camera pixel is co-registered with a scene pixel, so the resulting measurement value on each camera pixel is a linear combination of the spectral bands present at the scene pixel. The weights in this linear combination vary across the scene and are determined by a spatial filter placed in the imaging plane between the two dispersive 4-f lines, as explained in [17]. Dual disperser imagers are often used with reconfigurable spatial masks, for example digital micromirror devices (DMDs), so with different mask configurations, we make $N$ coded measurements of the scene.
Fig. 1. Simplified schematic of the dual disperser architecture and the programmable filtering $\textbf {H}$ implemented by the DMD mask.
For each region in the scene, a measurement $\textbf {m}_n$ corresponding to the values of the camera pixels in the region can be written as
For $N$ such acquisitions we can concatenate the results, giving
We obtain the panchromatic information $P_k$ via one acquisition with the spatial filter completely open. The spectrum $\textbf {s}$ for the homogeneous region can then be reconstructed using the measurements $\textbf {m}$ and this panchromatic information $P_k$ with a simple least squares procedure. However, as matrix $\textbf {G}$ can be ill-conditioned, we propose to use a Tikhonov regularization along the spectral dimension [18,19]. The estimated spectrum for the region is thus calculated by
From a computational point of view, the solution can be directly computed using Eq. (6), as the linear system to solve only has $W$ unknowns. Note that solution of the linear system has a constant computation cost regardless of the number of pixels $K$ in the region and the number of acquisitions $N$. Moreover, the computation costs of the matrix $\textbf {G}^{t} \mathbf {\Gamma }^{-1} \textbf {G}$ and the vector $\textbf {G}^{t}\mathbf {\Gamma }^{-1}\textbf {m}$ are linear with respect to $K$ and $N.$
Additionally, as each region can be computed independently, the algorithm is paralellizable, giving the potential for extremely fast reconstruction overall.
One important consideration for this algorithm is how the image can be divided into homogeneous regions, given no prior knowledge of the scene. We can consider that if two pixels have different spectra, there will likely be a difference in the panchromatic intensity, if we assume the presence of adjacent metamers is unlikely. For two adjacent but spectrally distinct regions, the difference in intensity is visible on the panchromatic as an edge, and thus some form of segmentation algorithm applied on the panchromatic image can be used to divide the image into regions. A gradual change in the panchromatic intensity does not pose a problem to the algorithm. Because the spectrum is constant over that region it can be reconstructed easily, although the segmentation algorithm may partition the area into separate regions. As segmentation is an ill-defined and non-trivial problem, we expect that this method will be most reliable for simple scenes that are easily segmented and obey the separability assumption e.g. a scene that is not highly textured, blurred and does not contain adjacent metamers or smooth spatial changes in spectra.
3. Experiment
The algorithm was tested experimentally using the dual disperser system described in [17], with an additional lens before the input. The scene has $W=110$ spectral bands in the range 425-650 nm, with 400 by 400 spatial pixels, and the spatial filtering is implemented using a DMD. Due to the orientation of the DMD, as well as size mismatch between the DMD mirror and pixel, the matrix $\mathbf {H}$ contains non-binary values. Whilst one would imagine that a fully 2D randomized DMD mask gives the best conditioning of Eq. (6), in practice the geometry mismatches lead to significant mixing of the data resulting in a less accurate reconstruction. To compensate for this effect, we are obliged to use a simpler mask design with only a 1D variation in the direction of dispersion. Consequently more acquisitions will be required than for an ideal 2D randomized mask, as there is less variation in the measurement data. The mask is designed according to [20], where long sequences of on or off mirrors are avoided to maximize variation in the information and improve conditioning of Eq. (6). This paper does not perform an extensive study into the ideal ratio of open mirrors (ROM), but we obtained good results with using $ROM=0.1$, with some scope for further optimization. The regularization parameter was chosen to be $\mu _\lambda =10^{3}$. A typical procedure to implement this algorithm consists of an acquisition of the panchromatic image, followed by one or more coded acquisitions with a 1D random mask.
3.1 Slit scanning
To evaluate the reconstruction method, ideally one would compare to the underlying ‘Ground Truth’ of the scene. However, for an experimental system this is not possible, so instead we compare to the reconstruction obtained via a straightforward slit scanning approach.
The data for both Slit Scanning (SS) and the Separability Assumption (SA) approaches are obtained with the same exposure time on the camera. The gray level on the DMD for the SS data was at its maximum value of 255. As SA combines several bands, the acquired images are brighter than for SS, so for the SA data the gray level of the DMD was reduced so as not to saturate the camera. We chose to keep the acquisition time constant to ease the comparison between SS and SA and avoid any nonlinearity at the camera level. In practice, the camera exposure time could be reduced for the SA data acquisition which would provide a further bonus with a reduction in terms of acquisition time. The SS reconstruction is limited by the dynamic range of the camera and its accuracy is limited for darker regions within the scene, where the signal to noise ratio is low, or where the signal is below the background level of the camera. Therefore, the SS data used for comparison to SA combines two data sets, with slit width of 1 mirror and slit width of 2 mirrors, to increase the dynamic range of the measurement. When the panchromatic intensity of the scene falls below a given threshold, the acquisition data from the wider slit is used to reconstruct the spectrum (further details in Supplement 1). The reconstruction with slit width of 2 mirrors will incidentally reduce the spectral resolution.
Similarly to compressed sensing approaches, the nature of the data acquisition also gives us the Fellgett’s advantage: compared to the SS acquisition, the SA acquisition has a higher SNR in darker regions.
To compare the two reconstructions we use the Spectral Angle Mapper (SAM), which measures the similarity of two spectra independent of their amplitude:
At a given position $(r,c)$ in the scene, the SAM is equal to zero when the spectra are identical, and equal to $\pi /2$ when the spectra are exactly orthogonal. The SAM is an ideal metric for this approach as it measures the spectral correspondence, rather than the spatial correspondence which is easily preserved by the geometry of the system.
3.2 Segmentation of the panchromatic image
Many segmentation algorithms are available, the latest of state of the art resorting to neural networks [21–24]. However, the focus of this paper is not on the segmentation algorithm itself, and we would ideally like that the method is robust enough that obtaining the exact segmentation is not crucial for the quality of the results (see Section 5). In practice, the panchromatic image is segmented into homogeneous closed regions using the classical watershed algorithm [25], implemented using MATLAB functions watershed. To alleviate noise, we first used an anisotropic diffusion algorithm [26], implemented using MATLAB function anisodiff2D, with gradient modulus threshold $\kappa =2\times 10^{-3}$, and 50 iterations. Anisotropic diffusion denoises the image, while preserving the contours, followed by a morphological gradient to improve the region-oriented segmentation, removing the gradients below a threshold $T,$ which varies the segmentation sensitivity. The segmentation algorithm returns the detected homegenous regions in the scene, which are bounded by a contour that is one pixel wide. The contours are likely to contain a spectral mixture of the two adjacent regions, but to focus primarily on the separability assumption they are discarded during the reconstruction of the HS image.
4. Results
4.1 Simple scene
We performed initial tests on a simple and easily segmented scene, consisting of a plastic brick wall illuminated with a white LED, shown in Fig. 2(a). The RGB image from the slit scanning (SS) reconstruction is shown in Fig. 2(b), and from the SA algorithm in Fig. 2(c), for $N=5$ coded acquisitions, the region contour highlighted in white. The DMD gray level was 255 for the SS acquisition and 140 for the SA acquisition. The panchromatic image was normalized between 0 and 1 and then segmented with a threshold of $T=1\times 10^{-2}$. Two examples of the reconstructed spectra are shown in Fig. 3, showing a good correspondence especially for the bright regions (Fig. 3(a)). Figure 4 displays the histogram and average of the pixel-by-pixel SAM values in the scene, dependent on the number of acquisitions. As one can see, only a few acquisitions are required for the results to converge to a good accuracy, with a low average SAM.
Fig. 2. (a) Panchromatic image, and RGB image for (b) slit scanning and (c) the SA method with $N=5$ acquisitions for the simple brick scene.
Fig. 3. Comparison of reconstructed SA spectra to the SS spectra, for pixels labelled (a) and (b) in Fig. 2(c).
The reconstructed data cube has a total of $Q =360$ regions, which means that the number of data points required to represent the scene is $(Q\times W)+(R\times C)=2 \times 10^{5}$, whereas the full datacube has size to $R \times C \times W=1.76\times 10^{7}$, pixels. This gives a total compression factor of 88 for the hyperspectral data cube. In terms of data acquisitions, the SA method requires $N+1=6$, compared to a normal slit scan with $W=110$ acquisitions (i.e. not high dynamic range), so the compression factor for the acquisition data is $W/(N+1)\simeq 18$. In terms of overall acquisition time for the whole sequence of images SA is 18 times faster. The reconstruction time for the SA method was 15 seconds, using non-optimized Matlab code.
This example also illustrates well the two main failure modes of the SA; firstly via regions that are too small to contain enough information, and secondly via adjacent metamers.
For very small regions (such as the bottom left of Fig. 2(c)), we can see on the RGB that the reconstructed spectra are inaccurate, as the reconstruction problem is under-determined, despite the regularization.
For adjacent metamers, Fig. 5 shows spectra for two different pixels in the same region, although by looking at the RGB image we can see that they are actually from different colored bricks (points (c) and (d) in Fig. 2(b)). When the scene contains two spectrally different zones, which are combined into a single region during segmentation, the reconstructed spectrum is a mixture of the spectra of the two individual zones. This is clear when we compare the spectrum constructed by SA (dotted red line in Fig. 5) to the average spectral value for that region found by SS (solid red line in Fig. 5).
Fig. 5. Comparison of reconstructed SA spectra to the SS spectra for pixels (c) and (d) labeled in Fig. 2(c), and for the average SS spectrum over the region.
4.2 Complex scene
Now we examine the performance with a more complex scene. The scene consists also of plastic bricks but has a 3D structure, gradients in light and dark due to shadows, and specular reflections, so is difficult to segment accurately. The DMD gray level was 255 for the SS acquisition and 125 for the SA acquisition. The RGB reconstruction of the scene using the slit scanning data is shown in Fig. 6(a). Figures 6(b), 6(c) and 6(d) show the reconstructed RGB results for the SA algorithm with varying segmentation thresholds.
Fig. 6. RGB images of the complex brick scene using (a) SS method, and (b-c) the SA method with $N=5$ acquisitions and segmentation threshold (b) $T=0.008$ (best segmentation in terms of average SAM), (c) $T=0.001$ (over-segmented), and (d) $T=0.02$ (under segmented).
Figure 7(a) shows the SAM histogram and average SAM as a function of the number of acquisitions, similar to the simple scene, although now as there are more small regions the convergence with $N$ is somewhat slower. For $N=5$, Fig. 7(b) shows that the reconstruction accuracy does not change much with the segmentation threshold, and there is a broad optimum region between $T=0.001$ and $T=0.015$. For an over segmented scene we would expect more small regions which do not contain enough information to be reconstructed accurately, whereas for the under segmented scene we are more likely to have two different regions merged into one, so the resulting spectra is a mixture of the two. Figure 7(c) shows how the average SAM of the region depends on the size of the region $K$, showing a large increase in the value and the variation in SAM as the region gets smaller. For a fully 2D mask we know the reconstruction problem can be over determined when $K>W/N=22$, i.e. for regions with more that 22 pixels we should have enough information to accurately reconstruct the spectrum. However, as we are using only a 1D mask, we should expect the limit on $K$ to be somewhat higher, depending on the shape of the region. Of course, regularization over the spectrum reduces the limit, depending on the smoothness of the spectra and the regularization parameter chosen. These results indicate that obtaining perfect scene segmentation is not vital to the success of the algorithm.
Fig. 7. Histogram of SAM values: (a) as a function or the number of acquisitions $N$ for $T=0.008$, and (b) as a function of the segmentation threshold for $N=5$; (c) For $N=5$, average SAM over a region as a function of the region size.
For the best segmentation, the reconstructed data cube has a total of 1138 regions, which means that the number of data points required to represent the scene is $(Q\times W)+(R\times C)=2.85 \times 10^{6}$, giving a total compression factor of 6.2 for the hyperspectral data cube. As for the simple scene, SA method permits a $\times 18$ reduction in the amount of data needed for reconstruction compared to the SS method. The reconstruction time for the SA method was 11 seconds (Matlab code).
5. Beyond panchromatic segmentation
As we have seen, the results are somewhat insensitive to the segmentation, and it may be possible to avoid panchromatic segmentation entirely. For example, Figs. 8(a) and 8(b) show the results for the SA approach when we arbitrarily divide the scene into $10\times 10$ or $20\times 20$ square regions. We can gain information on how well each region is reconstructed spectrally using the residual error, given by
Fig. 8. RGB reconstruction using SA methods on square regions or 'superpixels', region size (a) $10 \times 10$ pixels, and (b) $20 \times 20$ pixels, and $N=5$ acquisitions. (c) the average SAM plotted against the residual error $R$ of a region, for $10 \times 10$ regions.
6. Conclusion
The reconstruction algorithm is based on the particular properties of the dual disperser hyperspectral imager, namely the property of co-registration, and relies on spectral-spatial correlations in the scene. The algorithm reconstructs simple real-world scenes accurately and quickly using a small amount of data, also with an acquisition time much shorter than for a slit scanning acquisition. Full parallelization of the computation could theoretically reduce the computation time by a factor equal to the number of regions. Experimental improvements allowing a 2D mask, as well as investigation into the optimum ROM could reduce the number of acquisitions necessary or allow smaller regions to be reconstructed. Currently, due to the reliance on panchromatic segmentation, this method may not be well suited to challenging scenes that are difficult to segment e.g. highly textured scenes or with lots of metamers. However, preliminary tests using arbitrary division of the scene into regions show good results, and demonstrates that we can use the residual error as a metric on reconstruction accuracy. The use of the residual error gives scope to relevant areas of future research for this algorithm. The residual error could be used to avoid imperfect segmentation algorithms, or as an indicator for anomaly detection. Combined with the fast computation time, a dynamic and adaptive acquisition scheme could be developed, relying on the residual error or some other metric to inform subsequent acquisitions.
Funding
Agence Nationale de la Recherche (ANR-18-ASTR-0012-01).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. S. M. Nascimento, K. Amano, and D. H. Foster, “Spatial distributions of local illumination color in natural scenes,” Vision Res. 120, 39–44 (2016). [CrossRef]
2. A. Chakrabarti and T. Zickler, “Statistics of real-world hyperspectral images,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, (2011), pp. 193–200.
3. Y. Fu, Y. Zheng, I. Sato, and Y. Sato, “Exploiting spectral-spatial correlation for coded hyperspectral image restoration,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, (2016), pp. 3727–3736.
4. A. Wagadarikar, R. John, R. Willett, and D. Brady, “Single disperser design for coded aperture snapshot spectral imaging,” Appl. Opt. 47(10), B44–B51 (2008). [CrossRef]
5. A. A. Wagadarikar, N. P. Pitsianis, X. Sun, and D. J. Brady, “Video rate spectral imaging using a coded aperture snapshot spectral imager,” Opt. Express 17(8), 6368–6388 (2009). [CrossRef]
6. D. Kittle, K. Choi, A. Wagadarikar, and D. J. Brady, “Multiframe image estimation for coded aperture snapshot spectral imagers,” Appl. Opt. 49(36), 6824–6833 (2010). [CrossRef]
7. G. R. Arce, D. J. Brady, L. Carin, H. Arguello, and D. S. Kittle, “Compressive coded aperture spectral imaging: An introduction,” IEEE Signal Process. Mag. 31(1), 105–115 (2013). [CrossRef]
8. Z. Meng and X. Yuan, “Perception inspired deep neural networks for spectral snapshot compressive imaging,” in 2021 IEEE International Conference on Image Processing (ICIP), (2021), pp. 2813–2817.
9. L. Wang, T. Zhang, Y. Fu, and H. Huang, “Hyperreconnet: Joint coded aperture optimization and image reconstruction for compressive hyperspectral imaging,” IEEE Trans. on Image Process. 28(5), 2257–2270 (2018). [CrossRef]
10. H. Zhang, X. Ma, X. Zhao, and G. R. Arce, “Compressive hyperspectral image classification using a 3d coded convolutional neural network,” Opt. Express 29(21), 32875–32891 (2021). [CrossRef]
11. X. Yuan, T.-H. Tsai, R. Zhu, P. Llull, D. Brady, and L. Carin, “Compressive hyperspectral imaging with side information,” IEEE J. Sel. Top. Signal Process. 9(6), 964–976 (2015). [CrossRef]
12. L. Galvis, D. Lau, X. Ma, H. Arguello, and G. R. Arce, “Coded aperture design in compressive spectral imaging based on side information,” Appl. Opt. 56(22), 6332–6340 (2017). [CrossRef]
13. E. Salazar, A. Parada-Mayorga, and G. R. Arce, “Spectral zooming and resolution limits of spatial spectral compressive spectral imagers,” IEEE Trans. Comput. Imaging 5(2), 165–179 (2019). [CrossRef]
14. E. Salazar and G. R. Arce, “Coded aperture optimization in spatial spectral compressive spectral imagers,” IEEE Trans. Comput. Imaging 6, 764–777 (2020). [CrossRef]
15. P. Wang, J. Li, C. Qi, L. Wang, and J. Chen, “Input aperture restriction of the spatial spectral compressive spectral imager and a comprehensive solution for it,” Opt. Express 29(12), 17875–17889 (2021). [CrossRef]
16. M. E. Gehm, R. John, D. J. Brady, R. M. Willett, and T. J. Schulz, “Single-shot compressive spectral imaging with a dual-disperser architecture,” Opt. Express 15(21), 14013 (2007). [CrossRef]
17. E. Hemsley, S. Lacroix, H. Carfantan, and A. Monmayrant, “Calibration of programmable spectral imager with dual disperser architecture,” Opt. Commun. 468, 125767 (2020). [CrossRef]
18. J. Idier, Bayesian approach to inverse problems (John Wiley & Sons, 2013).
19. A. N. Tikhonov and V. Y. Arsenin, Solutions of ill-posed problems (V. H. Winston & Sons, 1977).
20. E. Hemsley, I. Ardi, S. Lacroix, H. Carfantan, and A. Monmayrant, “Optimized coded aperture for frugal hyperspectral image recovery using a dual-disperser system,” J. Opt. Soc. Am. A 37(12), 1916–1926 (2020). [CrossRef]
21. N. R. Pal and S. K. Pal, “A review on image segmentation techniques,” Pattern Recognition 26(9), 1277–1294 (1993). [CrossRef]
22. H. Narkhede, “Review of image segmentation techniques,” Int. J. Sci. Mod. Eng. 1, 54–61 (2013).
23. C. Amza, “A review on neural network-based image segmentation techniques,” De Montfort University, Mechanical and Manufacturing Engg., The Gateway Leicester, LE1 9BH, United Kingdom pp. 1–23 (2012).
24. N. Tajbakhsh, L. Jeyaseelan, Q. Li, J. N. Chiang, Z. Wu, and X. Ding, “Embracing imperfect datasets: A review of deep learning solutions for medical image segmentation,” Med. Image Anal. 63, 101693 (2020). [CrossRef]
25. J. B. T. M. Roerdink and A. Meijster, “The watershed transform: Definitions, algorithms and parallelization strategies,” Fundamenta Informaticae 41(1,2), 187–228 (2000). [CrossRef]
26. P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Trans. Pattern Anal. Machine Intell. 12(7), 629–639 (1990). [CrossRef]
