Probability of correct reconstruction in compressive spectral imaging

Coded Aperture Snapshot Spectral Imaging (CASSI) systems capture the 3-dimensional (3D) spatio-spectral information of a scene using a set of 2-dimensional (2D) random coded Focal Plane Array (FPA) measurements. A compressed sensing reconstruction algorithm is then used to recover the underlying spatio-spectral 3D data cube. The quality of the reconstructed spectral images depends exclusively on the CASSI sensing matrix, which is determined by the statistical structure of the coded apertures. The Restricted Isometry Property (RIP) of the CASSI sensing matrix is used to determine the probability of correct image reconstruction and provides guidelines for the minimum number of FPA measurement shots needed for image reconstruction. Further, the RIP can be used to determine the optimal structure of the coded projections in CASSI. This article describes the CASSI optical architecture and develops the RIP for the sensing matrix in this system. Simulations show the higher quality of spectral image reconstructions when the RIP property is satisfied. Simulations also illustrate the higher performance of the optimal structured projections in CASSI.


Introduction
Imaging spectroscopy requires sensing a large amount of spatial information across a multitude of wavelengths.The sensed signals are titled multispectral or hyperspectral images.Traditional imaging spectroscopy sensing techniques scan adjacent zones of the underlying spectral scene and merge the results to construct a spectral 3-Dimensional (3D) data cube.Push broom spectral imaging sensors, for instance, capture a spectral data cube with one FPA measurement per spatial line of the scene (Gehm et al., 2008).Spectrometers based on optical band-pass filters spectrally scan the source data cube (Gupta, 2008).These sensing techniques obey the well-known Nyquist criterion, which imposes a severe limit on the number of samples needed.More specifically, these methods require scanning a number of zones linearly in proportion to the desired spatial or spectral resolution.As the desired resolution increases, the required number of samples grows considerably, such that the cost of sensing a hyperspectral image is extremely high.Recently, with a mathematical technique for reconstructing sparse signals, called Compressive Sensing (CS), one can use fewer samples much lesser then predicted by the Nyquist criterion.This new technique involves diverse mathematical areas such as numerical optimization, signal processing, random matrix analysis, and statistics.The enormous potential of CS has been recently applied in areas such as microscopy, holography, tomography and spectroscopy (Arce et al., 2014;Brady et al., 2009;Studer et al., 2012;Wagadarikar et al., 2008;Yu and Wang, 2009).CS allows sensing a signal with a fewer number of samples than that required by the Nyquist criterion.Thus, CS has allowed overcoming a diverse number of physical sensing limits.This paper focuses on the application of CS in spectral imaging which is coined Compressive Spectral Imaging (CSI).CSI senses 2D coded random projections of the underlying scene, such that the number of required projections is far less than the linear scanning case.The Coded Aperture Snapshot Spectral Imaging (CASSI) system is an imaging architecture that effectively implements CSI.Thus, CASSI senses the 3D spectral information of a scene using 2D random projections.Figure 1   Figure 2 shows the main components of the CASSI architecture.In CASSI, the coding is applied to the image source density ƒ 0 (x, y, λ,) by the coded aperture T(x, y), where (x, y) are the spatial coordinates.The resulting coded density ƒ 1 (x, y, λ,) = T(x, y) ƒ 0 (x, y, λ,) is spectrally dispersed by a dispersive element before it impinges in the focal plane array (FPA), resulting in the coded field where T(x, y) is the transmission function representing ) is the optical impulse response of the system, and S(λ) is the dispersion induced by the prism.A linear dispersion by the prism is assumed.The compressive measurements across the FPA are obtained by the integration of the field ƒ 2 (x, y, λ ) over the detector spectral range.The discretized FPA measurements can be thus modeled as where L is the number of spectral bands sensed, is the discretized coded aperture with twodimensional structure, and ω jℓ is the noise system.
The single-shot CASSI system has been extended to a multi-shot system architecture such that multi-frame coded measurements are acquired as separate FPA measurements, each with a distinct coded-aperture pattern (Kittle et al., 2010).Equation ( 2) can be written in matrix form as y = Hf + ω, where y is the measurements vector, H is the matrix that represents the coded apertures and prism effect and f is the vectorized spectral data cube.Given the compressive measurements y, the objective of CS is to recover the signal f.Compressive Sensing (CS) dictates that an n-long signal f can be recovered from m ≪ n random projections y = H Ψ θ.Moreover, if it defines the m × n matrix A = HΨ as the sensing matrix, then the spectral data cube f = Ψθ is sparse in the base Ψ.The RIP in CASSI is thus fundamental to determine both the minimum number of shots needed for correct reconstruction and for the design and optimization of the coded aperture patterns (Wagadarikar, 2010).CS structured random matrices have been extensively analyzed for Toeplitz, random partial Fourier, and partial circulant matrices (Fornaiser, 2010).Recently a generalization of structurally random matrices has been developed (Do et al., 2012).The RIP in CASSI was first considered in Wagadarikar (2010), where it is assumed that the RIP in CASSI is satisfied.
The work in Wagadarikar (2010) then derives conditions on the coded apertures so that the RIP is better satisfied.These results, however, do not develop the explicit parameters for the bounds needed in the RIP, such as the probability of error, or the minimum number of FPA measurements.In this paper the RIP in CASSI is derived, the structure of the CASSI sensing matrix is formulated, and the RIP constants are expressed as a function of the structure of the random coded aperture patterns.A minimum number of random projections are needed so as to satisfy the RIP of the sensing matrix A (Candes and Tao, 2005).This paper establishes the theoretical probability of correct recovering in CSI.The results of this work can be used in the optimization of coded aperture patterns and to estimate a bound for the minimum number of shots needed for correct reconstruction.

CASSI Sensing Matrix
The goal in this section is to present the discrete model of the optical process in CASSI system in matrix form.Moreover, using this model allows carrying out the RIP analysis on the CASSI sensing matrix.The discretized output at the detector corresponding to the i th coded aperture T i j,ℓ for i = 1,!, K , where Y jℓ i is the intensity at the ( j,ℓ) position at the detector, whose dimensions are N × N + L − 1, F represents an N × N × L spectral data cube and ω jℓ is the white noise of the sensing system.In this work, the entries of the coded aperture are limited to T i j,ℓ varying terms in Equation ( 3), indexed by i are: first, the discrete coded apertures T i or its vector representation t i whose entries are realizations of a Bernoulli random variable with parameter p = P t ℓ i = { 1 } = 1− P t ℓ i = { 1 } and second, the number of shots K. (See Appendix A).
The i th FPA measurement Y jℓ i can be written in matrix notation as 3), where is the vector representation of the data cube F, and ω is the noise of the system.More specifically, the entries of f k can be expressed as where r = floor ℓ / N ( ) .The vectorization of the code aperture T i is defined as where r = floor ℓ / N ( ) .Similarly, the vectorization of the output Y i is defined as where r = floor ℓ / N ( ) .Using the above matrix representation, the output y i can be expressed as where diag (t i ) is an N 2 × N 2 diagonal matrix whose entries are the elements of the vectorized coded aperture t i , The matrix H i can be expressed as (Arguello and Arce, 2012) where 0 R is a N 2 × N(L − 1) zero-valued matrix, the structure of the V × V permutation matrix Θ v is The n × n matrix Θ n has the same structure of Θ v in Equation ( 10) and  Given the measurements y, the signal f is estimated by solving where τ is the regularization constant.Several algorithms have been developed to solve this l 2 − l 1 optimization.In this paper the Gradient Projection for Sparse Reconstruction (GPSR) algorithm is used (Figueiredo et al., 2007).

The RIP in CASSI
The sth restricted isometry property (RIP) δ s = δ s A ( ) of the matrix A ∈! KV ×n is defined as the smallest δ such that for all s-sparse vectors ϴ.Notice that equivalently Equation ( 13) is given by where the CASSI sensing matrix is A = HΨ, the matrix , and θ is a sparse representation of f in the base Ψ.Notice that A ∈! m×n where m = KV, and m ≤ n .Figure 3 and Figure 4 depict an example of the matrices H and A respectively for small values of N, L, and K.In practice, the dimensions of these matrices are much larger (Arguello and Arce, 2011).Figure 4 uses be rewritten as Given that the entries of the matrices A TT are Gaussian, then it is possible to apply the concentration inequalities developed by Ledoux to calculate the distribution of λ max in Equation ( 17) (Candes and Tao, 2005;Ledoux, 2001).More specifically, for each fixed t > 0 the concentration inequality is given by where σ max indicates the maximum singular value and o ( 1) is a small term tending to zero as K → ∞ which can be estimated as o 1 . Using Equation ( 25) and following a similar approach as that of Candes in (Candes and Tao, 2005), the RIP in CASSI can be written as where Pe = e Equation ( 17).
An important parameter of the CASSI system is the transmittance of the coded apertures.Algorithm 1 was developed in order to determine the relationship between the transmittance of the coded apertures and the RIP constant δ s .The objective of this section is to determine the transmittance that provides the best value of the constant δ s . The

Simulations
In order to verify the in Equation ( 26), an algorithm is proposed based on Monte Carlo Uniform Sampling.Also, two spectral data cubes are used to simulate the CASSI real system reconstructions.All experiments have been carried out in Matlab 2010b on a 3GHz Intel Core i7 and 32GB memory desktop computer.
Create CASSI matrix H i from Equation ( 9) and concatenate as Constant δ s vs the code aperture transmittance (parameter p), for various number of FPA measurement shots K, and S/n = 0,005.
The number of iterations in Algorithm 1 is defined as maxlter = n!/ (n − s) !s!.Since the number of all possible combinations of the sub-matrices is very large, the constant RIP δ s can be determined using Monte Carlo simulations.Therefore the index subspace S is sampled by using a uniform distribution (step 5) with size S ≤ s .The eigenvalues in all experiments are stored in the vector C. Finally, the probability of correct reconstruction is calculated from probability distribution of C by determining the number of eigenvalues smaller than 1.In this work, the parameters used in Algorithm 1 were N = 16, L = 8 , maxlter = 1 000 000 is the 2D-Wavelet Symmlet 8 basis and I ∈! L× L is the identity matrix.
pixels of spatial resolution and L = 33 spectral bands ranging from 400 nm to 720 nm.The third database is a color scene with 256 * 256 pixels of spatial resolution and L = 31 spectral bands ranging from 400 nm to 700 nm at 10 nm steps.In In Figure 10a, an example is shown where the spectral band selected is 500 nm.For the databases, the reconstruction is compared for different sparsity values using Bernoulli distributed coded apertures in CASSI with transmittance 50 %.
The reconstruction is realized using the GPSR algorithm (Figueiredo et al., 2007).Figure 6 shows the variation of δ s for different number of shots K, p, transmittances and S/n = 0,005.The constant δ s is minimum for P t ( ) ℓ = 1 } = p = 0,5 { , or equivalently when the transmittance of the code apertures is 50 %.Notice that a value of δ s = 1 indicates that reconstruction is not possible.
Similar performance is obtained for other values of S/n.The probability of reconstruction 1 − P e (see Equation ( 26)) is also compared to the results provided by a set of Monte-Carlo simulations.The results are shown in Figure 7, where a high probability depends on small sparsity and a higher number of shots.

Experiments with real data
In this section, the performance of CASSI sensing matrix is evaluated in the reconstruction of three spectral images (Foster et al., 2006).An example of the first database of 533 nm is shown in Figure 8a.This is a natural scene with 256 × 256 pixels of spatial resolution and L = 33 spectral bands ranging from 400 nm to 720 nm.An example of the second spectral database at 533 nm is shown in Figure 9a.The second spectral image is a city scene with 256 × 256

Minimum number of shots for reconstruction
While the Nyquist Shannon sampling theorem states that a certain minimum number of samples is required for sampling the spectral data cube, compressive sensing can reduce the number of measurements that need to be stored by the prior knowledge that spectral signal is sparse in a known basis and that sensing matrix satisfied the RIP. Figure 8d, Figure 9d and Figure 10d indicate the respective RIP constant δ s for correct reconstruction.The results indicate clearly that there are a minimum number of shots for correct reconstruction as a function of the sparsity for a noiseless case.The reconstruction of the band depicted in Figure 8a is shown in Figure 8b and Figure 8c for 4 and 9 shots respectively.Figure 8b illustrates the effect of using a number of shots less than that established by the RIP where the PSNR was 25 dBs with sparsity S/n = 0,065.For Figure 8c the PSNR was 31 dBs with sparsity S/n = 0,065.Notice that the number of shots used in Figure 8c is higher than the minimum recommended in Figure 8d.The number of shots used in Figure 8d is similar than the minimum recommended in Figure10d.The artifacts in the reconstructed image in Figure 8b, Figure 9b and Figure10c, with low PSNR, are clearly noticeable.On the contrary, Figures 8c, 9c and 10b show high quality in the reconstructed image.For Figure 10b and 10c the reconstruction performance are 28 and 24 dBs respectively, for the given number of shots 12 and 8.

Conclusions
The structured CASSI sensing matrix has been formulated and its RIP has been demonstrated.Simulations verified the theoretical RIP CASSI bounds.For optimal reconstruction, the coded aperture patterns must be designed with p = 0,5 or equivalent with a transmittance of fifty percent.Given a sparsity condition for signal representation, the minimum number of shots for reconstruction is determined by the variance of the representation basis and the variance of the random structure of the code apertures.An algorithm based on Monte Carlo simulations was developed to determine the optimal transmittance of the coded apertures.The experiments indicate the influence of the sparsity level and the number of shots in the performance of the system.
illustrates the CASSI optical architecture.The main components in CASSI are the coded aperture, the prism and the Focal Plane Array (FPA).The coded aperture patterns are the only varying elements in CASSI, the other optical elements remain fixed during the operation of the instrument.
times the usual matricial product, and 0 is a N 2 × N 2 (L − 1) zero valued matrix.

Figure 4 .
Figure 4. Structure of the matrix A for K = 2, N = 8 and L = 4 and Ψ is based on the 2D Symmlet Wavelet Transform.

Figure 3 .
Figure 3. Structure of the matrix H for K = 2, N = 6 and L = 3.The diagonals in the H matrix are the vectorized coded apertures.
r = S / n and H r ( ) = −r log r ( ) − 1− r ( ) log 1− r ( ) .Using Equation (26), the constant δ s is found as δ s = Estimate numerically for δ s with p and K To estimate numerically the restricted isometry constant δ s of order S from matrix A, Monte Carlo simulations are computed by sampling the sub-matrices A S T A s − I d ( ) from steps 1-3 in Algorithm 1 are used to construct the CASSI sensing matrix with parameters p and K.The coded apertures T i for i ∈ 0,…, K −1 } { are N × N symmetric matrices whose entries 1,−1 } { follow a Bernoulli distribution with probability p.The CASSI Sensing matrix A is formed by multiplying the CASSI matrix H with the sparsity basis Ψ.Subsequently, an iterative method is implemented to find the maximum eigenvalue of the possible sub-matrices A S T A s − I d ( ) .

Figure 5 .
Figure 5.The matrix depicts the matrix A TT for K = 2, N = 8, and S = 180.The off-diagonal entries of A TT can be seen as Gaussian entries.
The base representation Ψ diag Ψ D ,…,Ψ D ( ) where Ψ 2D is the 2D-Wavelet Symmlet 8 basis.The data cube F was approximate for different values of the sparsity S/n by removing the n − S less significant coefficients in magnitude from the sparse representation θ.The performance of matrices is analyzed under different sparsity levels.Figure 8d, Figure 9d and Figure10d indicate the performance (PSNR) analysis as a function of the sparsity and number of shots.

Figure 7 .
Figure 7. Probability of correct reconstruction 1 − P e vs the number of shots K and the sparsity S with transmittance p = 0,5.

Figure 8 .
Figure 8.(a) The slice of the spectral band 533nm.The respective reconstruction is illustrated for S/n = 0,065 for (b) 4 shots and (c) 9 shots.(d) PSNR as a function of the number of shots and sparsity.

Figure 9 .
Figure 9. (a) The slice of the spectral band 533 nm.The respective reconstruction is illustrated for S/n = 0,125 for (b) 12 shots and (c) 16 shots.(d) PSNR as a function of the number of shots and sparsity.

Figure 10 .
Figure 10.(a) The slice of the spectral band 500nm.The respective reconstruction is illustrated for S/n = 0,125 for (b) 12 shots and (c) 8 shots.(d) PSNR as a function of the number of shots and sparsity.