Multi-resolution reconstruction algorithm for compressive single pixel spectral imaging

Spectral imaging is useful in a wide range of applications for non-invasive detection and classification. However, the massive amount of involved data increases its processing and storing costs. In contrast, compressive spectral imaging (CSI) establishes that the three-dimensional data cube can be recovered from a small set of projections, that are generally captured in 2-dimensional detectors. Furthermore, the single-pixel camera (SPC) has been also employed for spectral imaging. Specifically, the SPC captures the spatial and spectral information in a single measurement. CSI reconstructions are traditionally obtained by solving a minimization problem using iterative algorithms. However, the computational load of these algorithms is high due to the dimensionality of the involved sensing matrices. In this paper, a multi-resolution (MR) reconstruction model is proposed such that the complexity of the inverse problem is reduced. In particular, this model uses the spectral correlation to group pixels with similar spectral characteristics. Simulation results show that the MR model improves the reconstruction PSNR in up to 9dB with respect to the traditional methods. In addition, the proposed approach is 79% faster, using only 25% of the measurements.


I. INTRODUCTION
Most imaging applications in areas such as astrophysics, environmental remote sensing, microscopy, surveillance, and biomedical image processing, often require high resolution images to discriminate between specific details of the scene.High-resolution sensors, however, increase the cost of the acquisition system.In particular, spectral imagery (SI), which consists on three-dimensional data sets with two spatial dimensions (x, y), and one dimension in the spectral domain λ, is traditionally employed in remote sensing for groundcover classification, mineral exploration, and agricultural assessment [1]; and biomedical imaging for noninvasive disease diagnosis and surgical guidance [2].The cost of acquiring SI using traditional scanning methods such as whiskbroom or pushbroom spectrometers depends on the desired resolution inasmuch as it determines the amount of scanned areas.In addition, SI captured with scanning-based systems involves massive amounts of data, which increases the cost of storing or processing.On the other hand, compressed sensing (CS) principles have been recently applied to SI acquisition [3], [4], this field has been called compressive spectral imaging (CSI).More specifically, CSI establishes that it is possible to retrieve a spectral image from a small number of samples, under the assumption that it has a sparse representation in some basis Ψ Ψ Ψ.In particular, a spectral image f ∈ R M •N •L has a dispersion level S, if it can be represented as a linear combination of S vectors on any basis Ψ, such that f = Ψθ with S M N L. Thus, instead of acquiring M N L samples, CS captures K M N L random projections of the scene, where K is not necessarily equal to the sparsity level S. The sensing process can be represented in matrix form as g = Hf , where H is the transfer matrix of the system [5].Because the number of measurements is considerably less than the number of voxels, the inverse problem given by f = H −1 g, is ill conditioned, leading to an infinite number of solutions.Therefore, reconstruction of the signal f is obtained using optimization algorithms that take advantage of the sparse representation of f in a transformation basis Ψ.Specifically, these algorithms obtain an approximation of θ and use the transformation f = Ψθ to get the desired spectral signal.The optimization problem is given by where τ is a regularization parameter.
In recent years, different optical architectures have been developed to implement the compressive sampling theory for the acquisition of spectral images, such as the spatial-spectral encoded compressive HS imager (SSCSI) [6], the coded aperture snapshot spectral imager (CASSI) [7], the dual-coded hyperspectral imager (DCSI) [8], the prism-mask multispectral video imaging system (PMVIS) [9], and the single pixel camera (SPC) [10].Each of these architectures presents its own advantages and drawbacks for different applications.However, most of them require 2D sensors.In contrast, SPC uses a single pixel detector, which provides reliable spectral images with lower cost hardware, since the single pixel detector can be a photodiode circuit or a point spectrometer (Whisk-broom).
As shown in Figure 1(a), this architecture is composed of the array lens L 1 and L 2 , a coded aperture T which is a blockunblock pattern implementable with a DMD, and a detector [11].
Traditional CSI reconstructions are obtained by solving Eq. 1 using iterative algorithms to retrieve a full-resolution approximation of the scene.In particular, each iteration of these algorithms calculates matrix products and inverses, however, the dimensionality of these matrices increases the computational load of the reconstruction process.In this paper, a MR reconstruction model is proposed such that the complexity of the inverse problem is reduced by employing superpixels, i.e. blocks of pixels.In particular, the blocks of pixels are  obtained by exploiting the spectral correlation of pixels with similar spectral characteristics.This model is based on the concept of multiresolution image (MRI) which assumes that blocks of pixels with similar characteristics can be grouped such that the number of elements in the full image is reduced [12], [13], [14].This approach has been previously applied for grayscale images with satisfactory results [13], [14], but usually the resolution for specific areas is selected by the user and does not take into account the actual characteristics of the image.It is worth noting that the MR approach has not been to date extended to CSI, the closest approximation to spectral imaging has been presented in [15], where the authors recover MR spectral images at the expense of high computational cost.

RECONSTRUCTION
The single pixel architecture illustrated in Figure 1(a) employs a coded aperture to modulate the input spectral information, then, a single point spectrometer (Whiskbroom) is used as a detector such that all the incoming modulated source is captured in a single measurement.Let a spectral image be F ∈ R M ×N ×L , where M and N are the spatial dimensions and L is the number of spectral bands.The sensing problem of a single spectral band f l can be modeled as g k l = h k f l , where l = 1, 2, . . ., L, h k is a row vector that contains all the physical phenomena behind the architecture including the coded aperture, and k indexes the captured snapshot.for each shot k.In general, the sensing model for all shots captured for the l-th band can be written as g l = Hf l , each of these employing a different coded aperture pattern, , H is the sensing matrix whose rows are the vectors h k .Furthermore, the sensing model for all the spectral bands and K shots is given by where g = (g 0 ) , Ĥ is the sensing matrix illustrated in Fig. 2 and can be obtained as a block diagonal matrix The compression ratio in this model is given by γ = K M N L , where the number of rows of Ĥ is equal to M N L. Notice that γ ∈ [0, 1].Full-resolution reconstructions from SPC measurements are obtained as in Eq. 1.To date, several optimization algorithms have been developed to solve the inverse CSI problem.In general, these algorithms work under the sparsity assumption of the underlying signal and their computational load is high due to the high dimensionality of the sensing matrices.In contrast, this work presents a reconstruction scheme that takes advantage of the spectral similarities of the pixels in the image.Furthermore, the scenes under analysis exhibit highly correlated areas in which several pixels can be grouped into blocks without losing information or inducing considerable errors.In particular, the proposed reconstruction scheme assumes that pixels of the same class have the same spectrum, such that they can be grouped into a superpixel, which reduces the amount of unknowns to recover in the inverse problem.Thus, the complexity of the resulting reconstruction problem is less than that of the full-resolution.Figure 1(b) summarizes the proposed reconstruction scheme, in which a low resolution (LR) reconstruction of the scene if first obtained; then, a high resolution version of the image is attained by interpolating the LR reconstruction which is then used as the input for a CS reconstruction algorithm in which the sensing matrix H is modified by a MR matrix D that accounts for the generation of the superpixels in the image.More specifically, in the MR reconstruction scheme, the sensing matrix H is composed by a Hadamard matrix W in conjunction with a decimation matrix D as in [16], and is given by where Z is a binary random matrix that helps to solve the inverse problem induced by the decimation matrix.The multishot sensing matrix, is then obtained as in Eq. 3. A fast LR reconstruction of f can be obtained by applying the inverse Hadamard transformation to the measurement set g.
Mathematically, the LR reconstruction can be obtained as ζ = ŴT g, where Ŵ is a block diagonal matrix, whose entries are W. Using the LR version of the scene ζ, a fast high-resolution (HR) approximation of f can be obtained by applying an upsampling operation to ζ without using an iterative algorithm.This HR approximation is later used to determine the superpixels of the image with respect to their spectra.Mathematically, the fast HR approximation of the scene is obtained as f = Uζ, where U is the employed upsampling operator, such as a bilinear interpolation matrix.
The correlation between the spectral signatures of the pixels in f determines the guidelines to group several pixels with similar spectral signatures into a superpixel, such that an estimation of the (MR) downsampling matrix D can be obtained.Specifically, the matrix D is a diagonal matrix whose entries correspond to the downsampling matrix for each shot ∆, then D is expressed as The procedure to generate the MR matrix is described in detail in Section III.Using D, H and the measurements g, the HR reconstruction problem can be rewritten as where τ is a regularization parameter and Ψ is the sparse representation basis.It is important to notice that ξ contains the values of the voxels in the HR image, and D indicates the spatial positions in which the elements of ξ should be placed.Thus, the reconstructed data cube f is given by f = DT ξ.Table I presents a summary of the proposed reconstruction methodology.

III. DESIGN OF MULTI-RESOLUTION MATRIX (∆)
The construction of the MR matrix is based on the analysis of the spectral signatures in the fast HR approximation of the scene as presented in Algorithm 1.The inputs for the algorithm are the fast HR approximation of Up-sampling with a bilinear interpolation (fast HR approximation) the scene f which is the only known information of the scene; the set of spatial coordinates of the image Ω Ω Ω = {(x, y) | x = 1, . . ., M − 1, y = 1, . . ., N − 1}; the spatial dimensions of the scene; the error tolerance σ; and the parameter for the block size η.In particular, the proposed method compares the spectral signature of each spatial point with respect to the average spectrum p of a block of pixels.The block size varies within the set T = 2 η , 2 η−1 , . . ., 1 .To start, the algorithm chooses a point (i, j) from Ω Ω Ω, and creates a block of size 2 η × 2 η and calculates the mean squared error (MSE) between the actual block B and a hypothetical block in which all points have spectrum p.If the maximum MSE is smaller than the fixed threshold σ then, the spatial points in the block are grouped into a single superpixel, thus the indicator vector Γ Γ Γ represents the pixels in the block and is given by where P is a replication of p whose size matches the size of B. In this case, the vector Γ Γ Γ is added as a new row of the MR matrix ∆ ∆ ∆.Furthermore, given that each point (x, y) can only belong to one superpixel, the set of pixels available for future blocks is given by Ω Ω Ω = Ω Ω Ω − B. On the other hand, if the maximum MSE is greater than the threshold, the algorithm randomly chooses another point (x, y) and performs the operations previously described.To keep track of the analyzed points, an additional set Ω Ω Ω is defined, such that after each evaluation, the pixel (x, y) is removed from the set by the operation Ω Ω Ω = Ω Ω Ω − (x, y).Then, before analyzing the MSE, the algorithm verifies whether the block can be constructed with the available pixels, mathematically, the new block should satisfy B = {(i, j)|(i, j) ∈ Ω Ω Ω}.Once all the pixels have been evaluated for the block size 2 η , i.e. | Ω Ω Ω |= 0, the algorithm moves to the next block size 2 η−1 , and runs the operations described above.The algorithm stops when all the pixels have been assigned to a block, i.e. | Ω Ω Ω |= 0. The size of the resulting matrix ∆ ∆ ∆ varies according to the data under analysis, however its number of rows is less than the original amount of unknowns to recover as shown in Eq.6.

V. CONCLUSION
A multi-resolution (MR) reconstruction scheme for spectral imaging recovery has been proposed.This approach exploits the spectral similarities within the image and reduces the complexity of the inverse problem.Simulation results show that the proposed method improves the reconstruction PSNR in up to 9dB and is 79% faster than the traditional reconstruction method, using a compression ratio of 25%.

Fig. 3 .
Fig. 3. Simulation results data cube 1 (a) RGB mapping original scene; (b) RGB mapping of the full-resolution reconstruction; (c) RGB mapping of the proposed reconstruction; (d) Blocks of pixels from proposed approach; (e) Recovered spectral bands.

Fig. 4 .Fig. 5 .
Fig. 4. Simulation results data cube 2 (a) RGB mapping original scene; (b) RGB mapping of the full-resolution reconstruction; (c) RGB mapping of the proposed reconstruction; (d) Blocks of pixels from proposed approach; (e) Recovered spectral bands.PSNR as a function of the number of recovered spectral bands for the three data cubes.The results are compared with the traditional full-resolution reconstruction.It can be noticed that the proposed approach improves the reconstruction quality in up to 9dB.In addition, Fig.7illustrates a comparison of the complexity of the reconstruction measured as the reconstruction time with respect to the number of recovered bands.It is worth noting that the measured time includes the generation of the MR matrix ∆ ∆ ∆, however this operation does not significantly increase the total time.These results show that

Fig. 6 .
Fig. 6.Comparison the PSNR of the reconstruction between proposed method compared with the traditional full-resolution reconstructions, is shown for (a)Target 1 (b)Target 2 and (c)Target 3.

Fig. 7 .
Fig. 7. Comparison of the time needly to reconstruct between proposed method compared with the traditional full-resolution reconstructions, is shown for (a)Target 1 (b)Target 2 and (c)Target 3.