1896

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

Colored Coded Aperture Design by Concentration of Measure in Compressive Spectral Imaging Henry Arguello, Member, IEEE, and Gonzalo R. Arce, Fellow, IEEE Abstract— Compressive spectral imaging (CSI) senses the spatio-spectral information of a scene by measuring 2D coded projections on a focal plane array. A 1 -norm-based optimization algorithm is then used to recover the underlying discretized spectral image. The coded aperture snapshot spectral imager (CASSI) is an architecture realizing CSI where the reconstruction image quality relies on the design of a 2D set of binary coded apertures which block-unblock the light from the scene. This paper extends the compressive capabilities of CASSI by replacing the traditional blocking-unblocking coded apertures by a set of colored coded apertures. The colored coded apertures are optimized such that the number of projections is minimized while the quality of reconstruction is maximized. The optimal design of the colored coded apertures aims to better satisfy the restricted isometry property in CASSI. The optimal designs are compared with random colored coded aperture patterns and with the traditional blocking-unblocking coded apertures. Extensive simulations show the improvement in reconstruction PSNR attained by the optimal colored coded apertures designs. Index Terms— Spectral imaging, coded aperture, CASSI, spectral selectivity, coded aperture optimization, rank minimization.

I. I NTRODUCTION

I

MAGING spectroscopy involves the sensing of a large amount of spatial information across a multitude of wavelengths. Traditional imaging spectroscopy sensing techniques scan adjacent zones of the underlying spectral scene and merge the results to construct a spectral data cube. Push-broom spectral imaging sensors, for instance, capture a spectral data cube with one FPA measurement per spatial line of the scene [1]. Spectrometers based on optical band-pass filters need to scan the scene by tuning the band-pass filters in steps [2]. The disadvantage of these techniques is that they require scanning a number of zones that grows linearly in proportion to the desired spatial or spectral resolution. In contrast, compressive spectral imaging (CSI) senses 2D coded projections of the underlying scene such that the number of required projections is far less than the linear scanning case. CSI exploits the fact that hyperspectral images can be Manuscript received June 28, 2013; revised December 17, 2013; accepted February 17, 2014. Date of publication March 5, 2014; date of current version March 18, 2014. This work was supported in part by ONR under Contract NOOO14- 1O-C-0-199 and in part by NSF under Grants EECS-0725422 and CCF-0915800. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Jose M. Bioucas-Dias. H. Arguello is with the Universidad Industrial de Santander, Bucaramanga 680001, Colombia (e-mail: [email protected]). G. R. Arce is with the Department of Electrical and Computer Engineering, University of Delaware, Newark, DE 19716-3130 USA (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TIP.2014.2310125

sparse in some basis representation. More formally, suppose a hyperspectral signal F ∈ R N×N×L , or its vector representation f ∈ R N·N·L , is S-sparse on some basis  , such that f =  θ can be approximated by a linear combination of S vectors from  with S  (N · N · L). Here, N × N represents the spatial dimensions and L is the spectral depth of the image cube. CSI allows f to be recovered from m random projections with high probability when m  S log(N · N · L)  (N · N · L). The Coded Aperture Snapshot Spectral Imager (CASSI) is one example of a CSI sensor that attains compressive measurements. The projections measured in CASSI are given by y = H f , where H is an N(N +L−1)×(N ·N ·L) matrix whose structure is determined by the coded aperture entries and the dispersive element effect. The CASSI architecture has been recently modified to admit multiple snapshots, each admitting a different coded aperture pattern. Multi-shot CASSI can thus yield a less ill-posed inverse problem and consequently improved signal recovering [3]–[6]. The time-varying coded apertures can be implemented using micro-piezo devices [3] or through the use of Digital Micromirror Devices [7]. Denote the t h FPA measurement as y = H f where H represents the effects of the t h coded aperture. The set of K FPA measurements, each with a different coded aperture is then assembled T  as y = (y0 )T , . . . , (y K −1 )T . The projections in CASSI can  θ = Aθθ where the matrix be alternatively expressed as y = H  is the sensing matrix. A = H The underlying data cube    θ ||2 + τ ||θθ ||1 is reconstructed as f˜ =  argmin||y − H θ  T where H = (H0 )T , . . . , (H K −1 )T , θ is an S-sparse representation of f on the basis  , and τ is a regularization constant [8]. Coded apertures in CASSI have been to date fabricated using materials such as chrome-on-quartz [5], [9] rendering coded aperture elements that are either opaque or transparent to the whole wavelengths of interest. These coded aperture masks are referred to as block-unblock coded apertures. New advances in microlithography and coating technology have allowed the design of patterned optical coatings which are used to create multi-patterned arrays of different optical filters. Indeed, coatings can be designed with high resolution patterns which are used on multispectral sensors [10], microelectro-mechanical devices [11], [12], optical waveguide-based devices, and gratings. Fig. 1 illustrates a typical colored coded aperture pattern. Colored coded apertures have been used recently in applications such as artificial compound eyes [13], multi-focusing and depth estimation [14], deblurring and matting [15].

1057-7149 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

ARGUELLO AND ARCE: COLORED CODED APERTURE DESIGN

1897

Fig. 2. Optical elements present in CASSI. The traditional block-unblock coded aperture is replaced by the colored coded aperture T (x, y, λ). Fig. 1. A colored coded aperture pattern is illustrated. Each pixel on the coded aperture is one of many possible optical filters whose spectral response can be selected.

This work incorporates the recent advances in patterned optical filters to the CASSI architecture such that traditional block-unblock coded apertures are replaced by multi-patterned arrays of selectable optical filters. The principal contribution of this paper lies on the optimal design of the parameters of the multi-patterned colored arrays. More specifically, the cutoff frequencies of the optical filters are optimized such that the Restricted Isometry Property of the CASSI is better satisfied. The optimal colored coded apertures reduce the number of 2D measurements required in CASSI for reconstruction. First, a detailed mathematical description of the structure of the matrix A is presented. The Hermitian matrix ATT AT is estimated where AT is a submatrix containing T arbitrary columns of the matrix A. It is shown that the entries of ATT AT can be modeled as sub-Gaussian random variables which are then used to determine the RIP constants in CASSI. Second, the sub-Gaussian constant is then obtained as a function of the colored coded aperture structure. The optimal colored patterns are calculated such that the sub-Gaussian constant is minimized and the respective distribution of the maximum eigenvalues of the matrices ATT AT is more concentrated [16], [17]. Third, two optimal colored coded aperture designs are presented. The first considers that the available colors are limited to low and high pass band filters. The second allows the incorporation of band pass filters in the design. A theoretical analysis compares the designed colored coded apertures with block-unblock coded aperture sets. Finally, simulations show the improvement attained by the colored coded apertures in terms of reconstruction Peak signal-to-noise ratio (PSNR). Also, the influence of the system noise in the reconstructed images is analyzed. II. CASSI W ITH C OLORED C ODED A PERTURES The CASSI system with colored coded aperture is depicted in Fig. 2 where the traditional block-unblock mask has been replaced by the colored coded aperture [5], [9]. The coding is realized by the coded aperture T (x, y, λ) as applied to the spatio-spectral density source f 0 (x, y, λ), resulting in the coded field f 1 (x, y, λ), where (x, y) are the spatial coordinates and λ is the wavelength. The coded density is then spectrally dispersed by a dispersive element before it impinges on the focal plane array as  T (x  , y  , λ) f0 (x  , y  , λ)h f 2 (x, y, λ) = 





×(x − S(λ) − x, y − y)d x d y



where T (x  , y  , λ) is the transmission function representing the coded aperture, h(x  − S(λ) − x, y  − y) is the optical impulse response of the system, and S(λ) is the dispersion induced by the prism. Each spectral slice of the data cube is thus spatially and spectrally modulated by the coded aperture and dispersed by the dispersive element [5], [9]. The compressive measurements across the FPA are realized by the integration of the field f2 (x, y, λ) over the detector’s spectral range sensitivity. The source f 0 (x, y, λ) can be written in discrete form as Fi j k where i and j index the spatial coordinates, and k determines the k t h spectral plane. The voxel Fi j k represents a measure of the intensity concentrated in a specific spatiospectral region of the source. Let  be the pixel pitch of the detector and let the discretization of the spectral axis be λk th for  k = 0, . . . , L − 1. The range of the k spectral band is λk , λk+1 where λk is the solution to the equation given by S(λk+1 ) − S(λk ) = ,

k = 0, . . . , L − 1.

(1)

f 0 (x, y, λ)d x d ydλ

(2)

The voxel Fi j k is then defined as λk+1 ( j+1) (i+1) 

Fi j k = λ

 k  =

j

i

f 0 (x, y, λ)d x d ydλ i j k

= ci j k · f 0 (x i , y j , γk ), for i, j = 0, . . . , N − 1, k = 0, . . . , L − 1 where, ci j k represents a quadrature weight, i j k represents the region defined by the intervals of integration, and (x i , y j , γk ) index values in i j k . Let Ti, j,k ∈ {0, 1} be the discretization of the colored coded aperture T  (x, y, λ) such that T (x, y,λ) =



Tij k rect

i, j,k



 x 1 1 y 1 λ −i − , − j − , − k − ,  2  2 σk 2

where it is assumed that the side length of the coded aperture square pixel is also equal to , and σk = λk+1 − λk . The horizontal and vertical spatial resolutions are limited by the resolution of the detector . The discretized coded aperture Tij k is a tridimensional structure whereas the traditional block-unblocking coded aperture is a two-dimensional structure. The t h discretized FPA measurement can be expressed as Yij =

L−1  k=0

Fi( j −k)k Ti( j −k)k + ωi j

(3)

1898

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

Fig. 3. Illustration of the spatio-spectral optical flow in CASSI. The q th slice of the data cube F is coded by a row of the coded aperture t and dispersed by the prism. The detector captures the intensity y by integrating the coded light. A colored coded aperture pixel labeled L, H, and B is colored with a low, high, or band pass filter respectively. A black pixel indicates a blocking pixel.

where Yij is the intensity at the (i, j )t h position of the detector whose dimensions are N × (N + L − 1). F is an N × N × L spectral data cube, and ωi j is the white noise of the sensing system. Fig. 3 depicts the physical sensing phenomenon in the colored CASSI for L = 6 spectral bands. There, the optical elements are represented by their effect on the discretized data cube. Each coded aperture pixel is colored with a low, high, or band pass filter. For instance, a green low pass element L means that the pixel permits the wavelengths λ0 , λ1 , and λ2 to pass through. A colored green high pass filter pixel H permits the wavelengths λ3 , λ4 , and λ5 to pass. The output Yij in (3) can be written in matrix notation as y = H f + ω

(4)

where y is a V -long vector representation of Yi, j with V = N(N + L − 1). f = vec([ f 0 , . . . , f L−1 ]) is the vector representation of the data cube F where f k is the vectorization of the k t h spectral band. The vectorization of F and Y , are detailed in Appendix B. The output y in (4) can be succinctly expressed as  ⎡ ⎢ ⎢ ⎢ ⎣

diag(t0 )



H

0 N(1)×N 2 diag(t1 )

 ⎤⎡ . . . 0 N(L−1)×N 2 ⎥⎢ ... ⎥⎢ ⎥⎢ .. ⎦⎣ .

0 N(L−1)×N 2 0 N(L−2)×N 2 . . . diag(tL−1 )

f0 f1 .. .

⎤ ⎥ ⎥ ⎥+ω ⎦

f L−1

(5) where diag(tk ) is an N 2 × N 2 diagonal matrix whose entries are the elements of the vectorized coded aperture planes tk detailed in Appendix B. If block-unblock coded apertures are used, then t0 = t1 =, . . . , = tL−1 . As will be shown shortly, the colored coded apertures lead to a richer structure of the matrices H inasmuch as their entries are less correlated. The ensemble of K measurements can be expressed as ⎡ ⎢ ⎢ y =⎢ ⎣

y0 y1 .. .

y K −1

⎤ ⎡ ⎥ ⎢ ⎥ ⎢ ⎥ =⎢ ⎦ ⎣

H

H0 H1 .. .

H K −1

⎤ ⎥ ⎥ ⎥ f + ω. ⎦

(6)

Fig. 4. The matrix H is presented for K = 2, N = 6, L = 3, and V = N (N + L − 1) = 48. White squares represent a one-valued element (unblocking light). Dotted line illustrates that the coded apertures are designed such that (tk0 ) j + (tk1 ) j = 1 for all j and k.

Multiple colored coded apertures can be implemented by coating a DMD or by moving a colored lithographic mask using a micro-piezo electric device. Fig. 4 depicts a random colored coded aperture-based matrix H for K = 2, N = 6, and L = 3. The white squares represent one-valued elements (unblocking light elements). As will be shown shortly, the coded apertures are designed such that all columns of H have the same number of white elements. Thus, (t0p ) j + (t1p ) j = 1 for all j and p. Notice in Eq. (5) and Fig. 4 that the i t h row of the matrix H contains up to L non-zero elements such that T  the i t h row of the matrix H = h0T , . . . , hTK V −1 can be written as ⎧   ⎨ t i k j i− V −k N , if i − i V = j − k j N and i j (hi ) j = ⎩0, otherwise, (7) for i = 0, . . . , K V − 1, and j = 0, . . . , N 2 L − 1, where k j =  j/N 2 , i = i /V , i ∈ {0, . . . , K − 1}, and N  = N 2 − N. In Eq. (7) and below, it is assumed that tkij =0 i−i V −k j N

when i − i V − k j N < 0. Remark that f is represented as f =  θ where θ is the sparse representation of the signal in the basis  . The sparsity

ARGUELLO AND ARCE: COLORED CODED APERTURE DESIGN

1899

basis can be expressed as  = W ⊗  2D where ⊗ is the Kronecker product. The basis function W is used to sparsify the wavelength axis, and the basis  2D is used to sparsify the spatial axes. Traditional basis functions are the wavelet transform, cosine transform, and pre-trained dictionaries  θ = Aθθ [4], [18]. Hence, y can be expressed as y = H f = H where A is referred to a sensing CASSI matrix which is a V K × N 2 L matrix. III. R ESTRICTED I SOMETRY P ROPERTY FOR C OLORED CASSI The Restricted Isometry Property of the matrix A is critical in the design of the coded apertures inasmuch as it guarantees the probability of correct reconstruction in CASSI. Definition 1 establishes the RIP in CASSI. Definition 1: Let the spectral data cube f be represented as f =  θ and assume |θθ | = S. The Restricted Isometry Property of the CASSI matrix A of order S is defined as the smallest δs such that (1 − δs ) ||θθ ||22

≤

||Aθθ ||22

≤

(1 + δs ) ||θθ ||22 ,

(8)

where the constant δs is given by δs =

max

T ⊂[ N 2 L ],|T |≤S

||A|TT | A|T | − I||22 ,

(9)

the operator ||.||22 is the squared norm from 2 into 2 , A|T | is a K V × |T | matrix whose columns are equal to |T | columns of A indexed by the set T , and I is an identity matrix [19], [20]. Defining the matrix A|T ||T | = A|TT | A|T | , Eq. (9) can be rewritten as   max λmax A|T ||T | − I (10) δs = T ⊂[ N 2 L ],|T |≤S where λmax (.) denotes the largest eigenvalue [20]. Let the entries of  be  j,k and let the columns of  be ψ 0 , . . . , ψ N 2 L−1 ], then using the structure of the matrices [ψ H in (7), the entries of A|T | can be written as   A|T | ir = hi ψ r L−1    = i−i V +k(N  ),r tk i =

i−i V −k N

k=0 L−1 



tki

m i −k N

k=0

m i +k(N  ),r

(11)

for i = 0, . . . , K V −1, r = 0, . . . , |T |−1, where i = i /V , i ∈ {0, . . . , K − 1}, m i = i − i V , i ∈ {0, . . . , V − 1}, 2 N  = N 2 − N, and r ∈ {0,  . . . , N L − 1}. Using (11), the entries A|T ||T | r,u can be expressed as K −1V −1  L−1  L−1   =0 i=0 k1 =0k2 =0

tk1



 i−k1 N

tk2



   Thus, E A|T ||T | r,r = 1 for r = 0, . . . , |T | − 1. The matrix A|T ||T | can be normalized by constraining the coded aperture  K −1  2 to satisfy = C, where C is a constant, =0 tk1 i−k1 N and by dividing the matrix A|T ||T | by C. This normalization procedure is explained in detail in Appendix A. Let the normalized matrix be A|T ||T | = A|T ||T | /C. Using Eq. (12) and defining the terms x,y,z =  K −1       ru 2D 2D ru =0 t x y t x−z y+z N , φ x,y =  x,r  x+y N,u , W x,y = Wx−y,u Wxr u [−y] + Wx,u Wx−yr u [y − 1] where u [x] is the step function of x, and defining the variable L−1−|q j |

β j ru =



k2 , p j ,q j Wkru2 ,q j /C,

(13)

k2 =0

  then the entries A|T ||T | can be succinctly expressed as r,u   A|T ||T |

r,u

=

NL 

φ ru p j ,q j β j ru

(14)

j =0

for r, u = 0, . . . , |T | − 1, where N L = N 2 (2L − 1) − 1, p j =  j/(2L − 1), and q j = j − p j (2L − 1) − (L − 1). A detailed procedure to obtain Eq. (14) is presented in Appendix C. It is assumed that the elements of the sparsifier matrices 2D  2D | < B and  2D and W are bounded such that |i,u 1 j,r |Wi,u W j,r | < B2 for all i, j, u and r . Thus, the variables φ ru such that p j ,q j and β j ru in (14) are bounded variables   ru  |φ p j ,q j | < B1 . Accordingly, each element A|T ||T | is r,u the sum of bounded random variables such that they can be modeled as sub-Gaussian random variables [20]. More  specifically, A|T ||T | ∼ Sub(ν 2 ) where the sub-Gaussian r,u constant ν is given by ν = max B1 r,u

NL 

β j ru .

(15)

j =0

The Restricted Isometry Property for sub-Gaussian random ensembles has been recently established in [21] and [22]. More specifically,   P ||A|T ||T | − I|| ≤ δs ≥ 1 −  (16) where  = 2(1 + 2/ρ) S e−δs (2−(1+ρ) ) K V c1 /ν with ρ = 2/ (e3 − 1) and c1 a constant independent of K and V . Notice that the probability of correct reconstruction 1 −  can be maximized by increasing the number of measurements K V or by minimizing the sub-Gaussian constant ν. 2

2 2

i+k1 N  ,ri+k2 N  ,u

i−k2 N

IV. C OLORED C ODED A PERTURE O PTIMIZATION (12)

for r, u = 0, . . . , |T | − 1. A necessary condition for the RIP of A is that the expected value of the diagonal elements of the submatrices A|T ||T | = A|TT | A|T | be equal to 1.

Equation (16) reveals that the RIP can be better satisfied by the sub Gaussian constant ν =  minimizing L maxr,u B1 Nj =0 β j ru , or equivalently by properly designing the coded aperture ensembles.

1900

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

Fig. 5. Illustration of an ensemble of 3 colored coded apertures. The number of bands is L = 5. A representative subregion of pixels in each shot is indicated with more intense colors which are taken into account in the optimization of the colored coded apertures. Each pixel is colored with a low, high, or band pass filter. The black pixels indicate blocking pixels.

A. Coded Aperture Optimization The following observations are taken into account in the optimization of the coded aperture ensembles: A.1) The sparsifying basis  2D and W are considered properly selected such that they induce high sparsity. Thus, the entries i,2Dj and Wi, j are considered fixed. Further, the unique assumption about the sparsifiers is that their entries are bounded, which is satisfied for most of the sparsifiers. Given that Wi, j and i,2Dj are assumed fixed, the optimization can  L be limited to minimize the variable ν  = Nj =0 β j r ∗ u ∗ , where u ∗ and r ∗ are the indices that maximize ν in (15). A.2) As illustrated in Fig. 3, each row vector of the coded aperture affects only one row at the detector, consequently each row of the coded aperture can be designed separately. The design procedure used for the i t h line can be repeated to design the other lines. simplification leads to optimizing  This  N−1  the variable ν  = 2L−2 k1 =0 v=0 β j r ∗ u ∗ , where j = (2L − 1) (i + v N) + k1 . By convenience i = 0 is chosen, such that j  = (2L − 1)Nv + k1 . A.3) It is only necessary to design an L-long section in each line of the coded aperture since the maximum spectral shifting produced by the dispersive element is L units. The other sections of the coded aperture   can be modulo repetition of the designed section such that tk j = tk j +L for all k and  [18]. This simplification is equivalent to the constrained optimization (minimization) of the variable 

ν =

+L−1 2L−2  v 0 k1 =0

β j r ∗ u ∗

(17)

v=v 0

where j  = (2L − 1)Nv + k1 and v 0 ∈ {L − 1, . . . , N (N − 1) − L + 1}. Fig. 5 illustrates an example of an ensemble of 3-colored coded apertures and L = 5 spectral bands. This figure illustrates the coded aperture sections taken into account to calculate (17) which are indicated with more intense colors. Each band is represented by the center wavelength λ j . B. Coded Aperture Optimization Principles The following observations can be used in the optimization of the coded apertures: B.1) The sum of correlated multiplications of the entries Wi, j leads to higher values of ν  in (17). More specifically,

Fig. 6. The spectral responses of the subregions of the coded apertures indicated in Fig. 5 are depicted. The gray squares indicate the spectral bands allowed through. Notice that the spectral responses of the filters are complementary among the shots such that the sum of the filters contains all the spectral components.

   L−1−|q j  | β j r ∗ u ∗ = k2 , p j  ,q j  {Wk2 −q j  ,u∗ Wk2 r∗ u −q j  + k2 =0  Wk2 ,u∗ Wk2 −q j  r∗ u q j  − 1 }/C in (17) tends to increase when k2 , p j  ,q j  = 0 for consecutive values of k2 . Thus, the value of ν  decreases as the randomness of the distribution of the non zero terms of k2 , p j  ,q j  increases. Let κ j = |L−1−q j  | k2 p j  q j  for j = 0, . . . , (2L − 1)L − 1, where k2 =0 j  = (2L−1)N j/(2L−1)+ j −(2L−1) j/(2L−1). Notice that κ j ∈ {0, . . . , L}. Further, let γ represent the frequency distribution of the values of κ , thus γi is given by γi =

(2L−1)L−1 

δ(κ j − i )/((2L − 1)L)

(18)

j =0

for i = 0, . . . , L, where δ(.) is the Kronecker delta function. The entry γi accounts for the probability that k2 , p ,q  = 0 for i consecutive values of k2 . Then, the optimization of the coded aperture aims at concentrating the distribution γi near to L i = 0. For instance, the minimization of the term i=1 γi ci leads to a concentration of the distribution of γL where ci are selectable positive cost weights such that i=1 ci = 1 and ci 1 < c for all i < i . A possible selection for ci is ci = i 1 2 2  L−1 α i where α is a selectable constant such that α > 1. α i / i=0 This selection of ci penalizes more the largest values of γi . B.2) Let the term U represent (17) when all the variables 2L−2 Wkru2 ,q j = 1. More specifically, the term U = k1 =0 v 0 +L−1  L−1−|q j  | k2 , p j  ,q j  /C. This term models the case v=v 0 k2 =0 when all the Wi, j are positive or all are negative (worst case scenario). It can be observed that the minimization of U tends to minimize ν  . Thus, the variable U is minimized instead of ν  . B.3) It can be observed that decreasing C increases considerably the value of U . Hence, the optimal value is C = 1. Fig. 6 depicts the spectral responses of the subregions of the coded apertures indicated in Fig. 5. The gray squares indicate the spectral bands that are let through by the respective filter. The constraint C = 1 is equivalent to designing the ensemble of the coded apertures such that the spectral responses of the filters are complementary among the shots in each spatial position of the ensemble. Accordingly, the sum of the spectral

ARGUELLO AND ARCE: COLORED CODED APERTURE DESIGN

1901

Fig. 9. Spectral responses of LH-Colored filter pixels and B-Colored filter pixels. (Left) three low pass and three high pass LH-Colored filter pixels are illustrated. (Right) three band pass colored filter pixels are shown.

Fig. 7. A geometric interpretation of optimization equation in (19). The columns of the matrices in Fig. 6 are circularly shifted according to the prism dispersion as indicated by the arrow. The term U is the total number of intersections of the gray squares between the columns of the shifted matrices. Three representative intersections are indicated with a green oval. The first, second, and third shifted matrix have 2, 8, and 3 intersections respectively such that U = 13.

Fig. 8. A shifted matrix of Fig. 7 is depicted. The intersection of the gray squares of the second and third column (dotted zone) of the indicated matrix increases the value of the term γ4 . This consecutive intersection pattern of the gray squares must be avoided in order to optimize the coded apertures.

responses through the shots spans the full spectrum under analysis. Fig. 7 illustrates a geometrical procedure to calculate the term U . The matrices in Fig. 6 are circularly shifted according to the dispersion effect of the prism. The sum of intersections of the gray squares between all the columns of each shifted matrix determines the value of U . Fig. 8 illustrates a poor design of the coded apertures. The consecutive intersection of the gray squares must be avoided in order to design properly the coded aperture ensemble. This consecutive intersection of points increases the value of γ4 . Given K the number of shots, the number of bands L, the data cube spatial dimension N,  the constant L−1 i i/ = α α > 1, the set of weights c i=0 α , the  K −1      i term x,y,z = =0 t x y t x−z y+z N , the term U = 2L−2 2L−1  L−1−|qk1 | k2 , pv ,qk1 , pv = Nv, qk1 = k1 − k1 =0 v=L k2 =0 (L − 1), and the function γi in (18), the optimization of the colored coded apertures can be succinctly expressed as argmin K −1 {t0 ,...,t L−1 }=0

subject to

L 

C. LH-Colored Coded Aperture Design The spectral response of low/high (L/H) pass colored coded apertures is constrained such that each coded aperture pixel can be a low or a high pass filter. The cut-off wavelengths of the filters are assumed to be selected from the subset {λ0 , . . . , λ L−1 }. Thus, there are 2λ L colored filters to be selected for each coded aperture pixel. More specifically, the spectral response of a λLj low pass colored coded aperture pixel is given by    1,  tk = i 0,

if k < λL j otherwise,

(20)

for k = 0, . . . , L − 1 where i ∈ 0, . . . , N 2 − 1, and  ∈ 0, . . . , K − 1. Similarly, the spectral response of a λH j high pass colored coded aperture pixel is given by    1,  tk = i 0,

if k ≥ λH j otherwise

(21)

for k = 0, . . . , L − 1. Notice that (20) can be alternatively written as     tk = δ k/λLj  (22) i

γ i c i + τ1 U

i=1 K −1  =0

response of the colored coded apertures can be constrained by cost and fabrication limitations. These limitations are given by the type of colored filters used in the fabrication of the coded apertures. This paper considers two specific designs. Firstly, the spectral response of the pixels in the coded aperture is limited to be either low or high pass filters. Secondly, band pass filters are permitted in the design of the coded apertures. The methodology used in these cases can be easily extended to other types of coded aperture limitations.

tk2

2 (L+k3 )N

=1

(19)

for k2 , k3 = 0, . . . , L −1, where τ1 is a regularization constant. The optimization in (19) can be seen as the minimization of the correlations of the spectral response of the ensemble of coded aperture pixels in a spatial range of length L. Further, the optimization aims at randomizing inevitable residual spectral correlations while maintaining constant the spectral energy of the ensemble at each spatial location. In practice, the spectral

where δ(·) is the Kronecker delta  function.  Equivalently, (21)   H can be written as tk i = 1 − δ k/λ j  . Fig. 9 (left) depicts the spectral response of three low pass and three high pass filters. Let the available set of low pass filters be L χ L = {λL 0 , . . . , λ L }, and the set of high pass filters H H χH = {λ0 , . . . , λ L }. The optimization in (19) is modified such that instead of designing the entries of the coded  } K −1 of cut-off wavelengths apertures, a set {ζ0 , . . . , ζ L−1 =0 is designed, where the cut-off wavelength ζ j ∈ {χ L ∪ χH} accounts for the spectral response along the axis k of the coded aperture pixel (tk ) L+ j N . More specifically, the optimization

1902

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

in (19) is rewritten as

V. S IMULATION AND R ESULTS L 

argmin  } K −1 {ζ0 ,...,ζ L−1 =0

γi ci + τ1 U,

i=1

subject to ⎧   ⎨δ k2 /ζ   ,   if ζk3 ∈ χ L  k3  tk2 = ⎩1 − δ k2 /ζ   , if ζ  ∈ χ H, (L+k3 )N k3 k3 K −1   =0

tk2

 (L+k3 )N

=1

(23)

for k2 , k3 = 0, . . . , L − 1. D. B-Colored Coded Aperture Design The spectral response of these colored coded apertures admits band pass filters. Thus, the spectral response of the coded aperture pixels can be selected from  a setHof band pass filters. Define the spectral response of a λL i , λi band-pass colored coded aperture pixel as    1,  tk = i 0,

H if λL i ≤ k ≤ λi otherwise

(24)

H for k = 0, . . . , L − 1 where λL i ≤ λi ∈ {0, . . . , L}. Again, Eq. (24)  can be written  usingHthe  delta Kronecker function /k δ k/λ  as tk i = δ λL i i . Fig. 9 (right) depicts the spectral response of three examples of band pass filters. The optimization in (19) is  rewritten such that a set of band pass filters is selected. Let ζk3 L , ζk3 H represent the cut-off wavelengths of the band-pass filter used to color the coded aperture pixel tk2 for k2 = 0, . . . , L − 1, then (19) (L+k3 )N is rewritten as L  argmin γi ci + τ1 U, K −1  ,ζ  ,...,ζ   {ζ0L 0H (L−1)L ,ζ(L−1)H }=0

i=1

subject to        tk2 = δ ζk3 L /k2  δ k2 /ζk3 H  K −1   =0

tk2



(L+k3 )N

(L+k3 )N

=1

(25)

for k2 , k3 = 0, . . . , L − 1. K −1  Evaluating each combination {ζ0 , . . . , ζ(L−1) }=0 in (23) requires about (2L − 1)L 2 (K + 1) multiplications. The total number of combinations in (23) is about (K 2L ) L such that the number of operations required search   for an exhaustive algorithm to solve (23) is O K 2L+1 L L+3 . Equivalently, the number of operations required  2 for an exhaustive search algorithm to solve (25) is O K L +2 L 3 . Given that the above optimizations involve only integer variables, and given that the number of required operations grows exponentially with L and K , a stochastic algorithm such as a Genetic Algorithm is used to solve (23) and (25) [4]. The GA reduces the complexity of the low-high pass filters optimization in (23)  to O M L K 2 where M is the number of selected individuals of the GA. Also, the optimization of the band pass filters in (25) is reduced to the order of O M K L 3 + 7M K ! .

In order to verify the CASSI colored coded aperture designs, a set of compressive measurements is simulated using the forward model in (3). These measurements are constructed employing a test spectral data cube which was acquired using a monochromator with wavelength steps of 1 nm. A bandpass filter (transmission window 450-670nm) is used to filter out unwanted spectral bands from the system. The image intensity was captured using a CCD camera AVT Marlin F033B, with 656 × 492 pixels, pixel pitch of 9.9μm, and 14 bits pixel depth. The resulting test data cube F has 256 × 256 pixels of spatial resolution and L = 16 spectral bands in the range 461nm to 596nm. The compressive sensing GPSR reconstruction algorithm [23] is used to recover the underlying data cube. The simulation of the model and reconstruction algorithms uses 64 bits arithmetic precision. The 2D spatial representation basis  2D used was the 2D-Wavelet Symmlet 8 basis. The spectral axis sparsifier W was the 1-D cosine transform basis which has the advantage of low computational cost compared to 3D wavelet transforms. The simulations are performed in a desktop architecture with an Intel Core i7 3.3GHz processor, 32GB RAM, and using Matlab R2012. A. Code Design Results A Genetic Algorithm (GA) was used to solve the optimization problem in (23) and (25) for K = 1, . . . , 14, L = 16, and N = 256. The parameters used for the GA are: number of individuals 36, number of generations 400, mutation rate 0.004, and crossover rate 0.3. The GA takes in average 2 hours to find an optimal solution and converges 100% of the time. Fig. 10 shows a realization of the optimal LH-Colored coded apertures for K = 3 shots. There, a 16 × 16 spatial section of the optimal coded apertures is presented. The cut-off wavelength of the filter associated with each pixel of the coded aperture in Fig. 10 is represented by the respective color. A small black square in the center of a pixel indicates that that filter is high pass. A black pixel indicates that it does not let pass the light. Notice in Fig. 10 that the spectral response of the coded apertures along several shots is complementary, such that the entire set of frequencies is sensed in each spatial position. Fig. 11 depicts the distribution of the cutoff wavelengths across pixels. Notice in Fig. 11 that the distribution is centered around 3 points, the central wavelength (L/2), the 0 wavelength (block) and the L wavelength (unblock). Fig. 11 suggests that a coded aperture ensemble with block, unblock, and LH-colored coded apertures centered at the middle wavelength can be the optimal selection when the available set of cut-off wavelengths is limited. The optimized Band-Pass colored coded apertures are illustrated in Fig. 12 for K = 3 shots. In this case, each optical filter pixel is represented by a square, where the colors indicate the wavelengths in the band pass of the pixel. An example of the spectral response of two colored pixels is presented at the bottom of Fig. 12. Fig. 13 illustrates the distribution across pixels of the band pass cut-off wavelengths of the

ARGUELLO AND ARCE: COLORED CODED APERTURE DESIGN

1903

Fig. 12. Geometric interpretation of colored coded apertures for B-filters (3 shots) is illustrated. An 8×8 pixel portion of the coded apertures is shown. The colors inside each pixel indicate the wavelengths that are permitted to pass in the respective pixel. The spectral response of two representative pixels is also shown.

Fig. 10. A geometric interpretation of the colored coded apertures for LH-Colored filters (3 shots). A small black rectangle in the middle of a pixel indicates a high pass filter pixel. If the small rectangle is not present, the pixel represents a low pass filter. In both cases, the color of the pixel represents the cut-off wavelength. A black pixel indicates a blocking pixel. An example of the wavelength response of two pixels is also depicted.

Fig. 13.

  Cut-off wavelength λL , λH pair occurrence probability for the

coded aperture filters shown in Fig. 12.

Fig. 11. Distribution of cut-off wavelengths for (a) Low pass filters, and (b) High pass filters of the coded apertures in Fig. 10.

designed coded apertures. The y-axis indicates the lowest cutoff wavelengths λL and the x-axis indicates the highest cut-off wavelengths λH of the band pass filters. Fig. 13 suggests that the most probable filters are the low band pass filter (cutoff L/3), high band pass filter (cut-off 2L/3), and band pass filters with cut-off wavelengths (L/3, 2L/3). When the cut-off wavelengths of the band pass filters are limited, the indicated filters are likely to be the optimal. B. Reconstruction from LH-Colored Coded Aperture Measurements The coded aperture ensembles designed in the previous section are used in capturing the CASSI compressive measurements and in the corresponding image reconstructions.

Fig. 14. Mean PSNR of the 256 × 256 × 16 reconstructed data cubes using Block-Unblock, Random LH-Colored and LH-Colored coded aperture ensembles.

Fig. 14 depicts the PSNR of the reconstructed images as a function of the number of shots for the optimized LH-Colored coded aperture ensembles. The results obtained with an ensemble of coded apertures based on a random selection of the cut-off wavelengths of the filters (LH-Random Colored) is also presented in Fig. 14. For comparison purposes,

1904

Fig. 15. The resulting spectral data cubes are shown as they would be viewed by a Stingray F-033C CCD Color Camera. (a) The original data cube and its reconstruction using (b) Optimal LH-Colored coded apertures, (c) Random LH-Colored coded apertures, and (d) Block-Unblock coded apertures.

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

Fig. 16. Zoomed versions of the images in Fig. 15: (a) Original data cube, and reconstructions by (b) Optimal LH-Colored coded aperture, (c) Random LH-Colored coded aperture, and (d) Block-Unblock coded apertures.

the reconstruction PSNR for a random block-unblock coded aperture ensemble is also included in Fig. 14. In the latter, the pixels are realizations of Bernoulli random variables with p = 0.5. To better visualize the reconstructed data cubes. The whole reconstructed data cube is presented in Fig. 15 as it would be viewed by a Stingray F-033C CCD Color Camera. There, it is possible to observe the better performance of the optimized LH-Colored coded aperture ensembles where a PSNR improvement of up to 4.56 dB is attained for 3 shots compared with the LH-Random Colored coded aperture ensembles. Fig. 16 shows a zoomed version of Fig. 15. The quality of reconstruction along the spectral axis is compared for the different coded aperture ensembles in Fig. 17 where three representative spatial pixels are considered. Fig. 17 indicates the higher reconstruction performance of the optimized LH-Colored design over the random structures. C. Reconstruction from B-Colored Coded Apertures Fig. 18 depicts the mean PSNR of the reconstructed images when the optimized B-Colored coded aperture ensembles are used. The results for a random selection of the band pass wavelengths and the block-unblock coded aperture ensembles are also included in Fig. 18. A zoomed version of a selected region of the original 6t h spectral band with center wavelength 506nm, is presented in Fig. 19. There, the reconstructions for 3 shots using the different coded aperture ensembles are shown. Additionally, Fig. 20 illustrates the whole reconstructed data cube as it would be viewed by a Stingray F-033C CCD Color Camera. A zoomed version of this figure is depicted in Fig. 21, where it can be observed that the optimized

Fig. 17. Reconstruction along the spectral axis of three spatial pixel locations as indicated in (a). The spectral responses are illustrated for (b) the pixel B, (c) the pixel C, and (d) the pixel D.

B-Colored coded aperture ensemble provides an improvement of up to 9dB in PSNR over the random selection of the band pass filter wavelengths. Furthermore, the PSNR reached by the B-Colored design is 4.6 dB higher than the LH-Colored design for the illustrated three shots reconstructions. The quality along the spectral axis is depicted in Fig. 22 where three representative spatial pixels are considered. Fig. 22 indicates the higher spectral reconstruction

ARGUELLO AND ARCE: COLORED CODED APERTURE DESIGN

1905

Fig. 18. Mean PSNR of the reconstructed data cubes using Block-Unblock, Random B-Colored and B-Colored coded apertures ensembles.

Fig. 20. The resulting spectral data cubes are shown as they would be viewed by a Stingray F-033C CCD Color Camera. (a) The original data cube, and its reconstruction by using (b) Optimal B-Colored coded apertures, (c) Random B-Colored coded apertures, and (d) Block-Unblock coded apertures.

Fig. 19. Zoomed versions of a selected region of the 6th spectral band. (a) Original, and reconstructions using (b) Optimal B-Colored coded aperture (45.5 dB), (c) Random B-Colored coded aperture (40 dB), and (d) BlockUnblock coded apertures (35.9 dB).

performance of the optimized B-Colored design over the random colored coded aperture structures. D. Noise Analysis and General Results Several experiments were conducted to verify the stability of the reconstructions. The consistence in the reconstructions with respect to the GA, the effect of sparsity of the scene, noise of the measurements, extinction ratios of the optical filters, and misalignment of the pixels are studied. Gaussian noise with zero mean was added to the measurements. The Signal to Noise Ratio (SNR) is calculated according to S N R =  10 log10 μy /σnoise , where μy is the mean value of the measurement y and σnoise is the standard deviation of the signal noise. This equation is used as the signals analyzed are digital images represented by positive pixel values. D.1) Consistency of Reconstructions: The first set of experiments considers the consistency of the genetic algorithm solutions. Since the GA is a stochastic optimization method, different realizations of the GA generate distinct coded

Fig. 21. Zoomed versions of the images in Fig. 20: (a) Original data cube. Reconstructions by (b) optimal B-Colored coded aperture, (c) Random B-Colored coded aperture, and (d) Block-Unblock coded apertures.

aperture ensembles. Ten experiments were performed for each number of shots and for each kind of filter ensemble. Table I indicates the maximum, minimum, and median PSNR of the reconstructed images for the different coded aperture structures. The PSNR values obtained with B-colored coded apertures are higher than those obtained with random band

1906

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

Fig. 23. Mean PSNR of the reconstructed data cubes for different values of SNR in the measurements (a) 15dB, (b) 20dB, (c) 30dB, (d) 40dB.

Fig. 22. Reconstruction along the spectral axis of three spatial pixel locations as indicated in (a). The spectral responses are illustrated for (b) the pixel B, (c) the pixel C, and (d) the pixel D.

TABLE I M EAN R ECONSTRUCTION PSNR IN dB

Fig. 24. The phase transition diagrams are illustrates for (a) suboptimal design based on Band-pass colored coded apertures, and (b) optimal B-Colored coded apertures.

pass filters (Random B-Colored). This relationship also holds between LH-colored and random Low-Pass and High-Pass colored coded apertures. The lowest values in PSNR are obtained with the Block-Unblock coded apertures. D.2) Noise in Measurements: Fig. 23 shows the mean PSNR of the reconstruction for the different coded aperture ensembles as a function of the number of shots for different SNR levels in measurements. Again, ten experiments were performed for each type of colored coded aperture. Despite the loss in reconstruction quality, the optimized ensembles conserve better performance than the random structures. D.3) Sparsity Analysis: In order to determine the effect of the sparsity in the random and optimal designs phase transition diagrams are shown in Fig. 24(a) and (b) respectively. The optimal design presents a higher performance for the different sparsity ratios S/N 2 L where S is the number of non-zero elements of the underlying data cubes.

D.4) Impact of the Extinction Ratios of the Optical Filters: The optimization of the colored coded aperture considers that the filter are transmissive in a given zone and completely opaque elsewhere (ideal). In practice, the optical filters are not completely opaque in the blocking region. Fig. 25(a) illustrates the impact of the extinction ratios of the optical filters. Fig. 25(a) shows the performance for two commercially available spectral responses of the optical filters named order1 and order2 respectively. It can be observed that extinction ratios affect the performance of the optical filters, however the gain of the optimal colored coded aperture remains over the traditional CASSI. D.5) Coded Aperture Misalignment Analysis: Simulations considering the misalignment between the features of the coded aperture and the pixels of the detector are illustrated in Fig. 25(b). The misalignment is expressed as a percentage of the size of the feature pixel of the detector. Notice that Fig. 25(b) considers simultaneously the effect of the extinction ratios and the misalignment. It can be observed that the colored coded apertures preserve the gain over the traditional blockunblock CASSI systems. D.6) Singular Value Spectrum Analysis: The singular value spectrums for the sensing matrices based on the random selection of the filters (Random LH and random B-colored), the optimal designs (LH colored and B-colored), and the traditional block-unblock CASSI systems are presented in Fig. 26.

ARGUELLO AND ARCE: COLORED CODED APERTURE DESIGN

Fig. 25. (a) The impact of the extinction ratios in the optimization of the colored coded apertures is illustrated. Three different spectral responses are simulated (Ideal, order1 and order2) for the optimal (B-Colored) and the random designs (Random B-Colored). The traditional block-unblock CASSI is also indicated. (b) The impact of the pixels misalignment is presented for three different percentages of misalignment. Each percentage of misalignment is simulated for the three extinction ratios indicated in (a).

1907

the improvements range from 1 dB to 23 dB. A similar behavior is experienced by the optimized band pass structures, where the improvements range from 1dB to 26 dB. The results indicate that the optimized structures lead to a better reconstruction along the spectral axis. The optimized coded aperture ensembles also exhibit a better performance even when the measurements are contaminated with additive noise, the optical filters present extinction ratios, and the pixels of the coded aperture have misalignment. Ongoing work aims at the optimization of colored coded apertures when the number of colors is much lower than the number of spectral bands. A PPENDIX A N ORMALIZATION OF THE M ATRIX A|T ||T | The diagonal elements of the matrix A|T ||T | in (28) can be expressed as 

A|T ||T |

where γ =

Fig. 26. The singular value spectra for the sensing matrices based on the random selection of the colored filters (Random LH and random B-colored), the optimal designs (LH-colored and B-colored), and the traditional random block-unblock are illustrated. The condition number of the matrices λ1 /λr is also indicated.



=

r,r

−1  L−1 K −1V   =0 i=0 k1 =0

tk1

2

2 i+k + γ (26)  1 N ,r

i−k1 N

K −1V −1    tk1 t i+k1 N ,ri+k2 N ,r . i−k1 N k2 i−k2 N =0 i=0 k1 =k2

Let C be a selectable constant and  let the ensemble of coded K −1  2 apertures be constrained to satisfy =0 tk1 i−k N = C, for 1 all i, k1 , then (26) can be written as ⎞ ⎛   V −1 bi A|T ||T | r,r 1 ⎝  2 = Ci+k + γ⎠  1 N ,r C C i=0 k1 =ai

There can be observed that all random-type designs exhibit a similar behavior. The optimal designs, on the other hand, present a different distribution of the singular values. It can be observed that the condition number of the optimal designs λ1 /λr is significantly smaller than that of the random designs. Thus, the optimal design leads to well condition sensing matrices. Notice also that the matrices based on the optimal designs do not admit truncated SVD approximations inasmuch as these matrices are almost full rank. VI. C ONCLUSION The mathematical model of the Colored Coded Aperture Compressive Spectral Imaging system is presented. The restricted isometry property in colored CASSI and therefore the probability of correct reconstruction is maximized by properly designing the spectral response of optical filters representing the colored coded aperture pixels. The optimization is formulated for two scenarios, the first constraints the optical filters to be low and high pass filters, and the second scenario admits band pass filters. The random low and high structures are superior than the block-unblock coded apertures. The optimized low and high pass structures are superior than the random colored patterns. For the optimized colored coded apertures, the PSNR increases rapidly with the number of measurement shots. For instance, in this case



A|T ||T |

 r,r

=

2 −1 L−1 N 

i=0 k1 =0

= 1+

2 i+k

1N

2 , r

1 γ C

+

1 γ C (27)

 2 where the condition tk1 i−k N = 0 when i − k1 N < 0, it is 1 taken into account by the summation limits ai and bi where ai = 0, bi = L if |i | ≥ L, ai = 0, bi = mod N i if i < L, and ai = mod N |i | + 1, bi = L if −i < L. Equation (27) uses the fact that the sparsity basis function is orthonormal, thus  N 2 L−1 2  N 2 −1  L−1 2  j,r = 1. The term γ i=0 k1 =0 i+k1 N 2 ,r = j =0 in (27) has the same structure of the off-diagonal elements of the matrix A|T ||T | such that all the elements of the matrix A|T ||T | − I have the same statistical structure. A PPENDIX B V ECTORIZATION OF THE C OLORED CASSI VARIABLES The entries of f k can be expressed as ( f k ) j = F( j −r N)rk , for j = 0, . . . , N 2 − 1, k = 0, . . . , L − 1  where r = Nj . Let the vectorization of y in (4) be defined as   y j = Y(j −r N)r , for j = 0, . . . , V − 1,  = 0, . . . , K − 1

1908

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

 where r = Nj and V = N(N + L − 1). Let the vector tk represent the vectorization of the 3D coded aperture Tij k at the 2D k t h plane. Thus, the elements of tk are given by   for j = 0, . . . , N 2 − 1 tk j = T(j −r N)rk ,  where r = Nj . A PPENDIX C C ALCULATION OF THE E NTRIES OF THE M ATRIX A|T ||T |   Using (11), the entries A|T ||T | r,u can be expressed as −1  L−1  L−1 K −1V  

tk1

=0 i=0 k1 =0k2 =0



 i−k1 N

tk2



i+k1 N  ,ri+k2 N  ,u

i−k2 N

(28) for r, u = 0, . . . , |T | − 1. To simplify the notation, when i − k1 N < 0 or i − k2 N < 0 in (28), it is again assumed that tk1 i−k N = 0 and tk2 i−k N = 0. Given that the sparsity 1  2 basis is  = W ⊗  2D , A|T ||T | r,u in Eq. (28) can be expressed as L−1  L−1 K −1V −1    =0 i=0 k1 =0k2 =0

2D 2D βik1 k2 i−k  Wk1 r i−k N, Wk2 ,u 1 N, 2 r



tk1



u

(29)

  tk2

, x = x i−k1 N i−k2 N x x  N 2 , and x =   . Define the terms x,y,z N2 NK2 −1       ru 2D 2D ru =0 t x y t x−z y+z N , φ x,y =  x,r  x+y N,u , and W x,y Wx−y,u Wxr u [−y] + Wx,u Wx−yr u [y − 1] where u [x] where βik1 k2

=

−

= = is the step function of x. Reorganizing the terms in (29) and taking into   account the mentioned normalization, the entries A|T ||T | can be written as r,u

  A|T ||T |

r,u

L−1 

=

N−1  N−1 

k1 =−L+1 i=0 j =0

×

L−1−|k  1| k2 =0

φ ru j N+i,k1

k2 , j N+i,k1 Wkru2 ,k1 /C.

(30)

Let N L = N 2 (2L − 1) − 1, p j =  j/(2L − 1), and q j = j − p j (2L − 1) − (L − 1), then (30) can be written as L−1−|q j | NL      k2 , p j ,q j Wkru . A |T ||T | r,u = φ ru p j ,q j 2 ,q j j =0

(31)

k2 =0

Define the variable β j ru , depending on the coded aperture ensemble, as L−1−|q j |

β j ru =



k2 , p j ,q j Wkru2 ,q j /C,

(32)

k2 =0

  then the entries A|T ||T |

r,u



A|T ||T |



for r, u = 0, . . . , |T | − 1.

r,u

can be expressed succinctly as =

NL 

φ ru p j ,q j β j ru

j =0

(33)

R EFERENCES [1] R. G. Sellar and G. D. Boreman, “Classification of imaging spectrometers for remote sensing applications,” Opt. Eng., vol. 44, no. 1, pp. 013602-1–013602-2, 2005. [2] G. Nahum, “Imaging spectroscopy using tunable filters: A review,” Proc. SPIE, vol. 4056, pp. 50–64, Apr. 2000. [3] D. Kittle, K. Choi, A. A. Wagadarikar, and D. J. Brady, “Multiframe image estimation for coded aperture snapshot spectral imagers,” Appl. Opt., vol. 49, no. 36, pp. 6824–6833, Dec. 2010. [4] H. Arguello and G. R. Arce, “Code aperture optimization for spectrally agile compressive imaging,” J. Opt. Soc. Amer. A, vol. 28, no. 11, pp. 2400–2413, Nov. 2011. [5] H. Arguello, H. Rueda, Y. Wu, D. W. Prather, and G. R. Arce, “Higherorder computational model for coded aperture spectral imaging,” Appl. Opt., vol. 52, no. 10, pp. D12–D21, 2013. [6] 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., vol. 31, no. 1, pp. 105–115, Jan. 2014. [7] Y. Wu, I. O. Mirza, G. R. Arce, and D. Prather, “Development of a digital-micromirror-device-based multishot snapshot spectral imaging system,” Opt. Lett, vol. 36, no. 14, pp. 2692–2694, 2011. [8] H. Arguello, C. Correa, and G. R. Arce, “Fast lapped block reconstructions in compressive spectral imaging,” Appl. Opt., vol. 52, no. 10, pp. D32–D45, 2013. [9] A. A. Wagadarikar, R. John, R. Willett, and D. Brady, “Single disperser design for coded aperture snapshot spectral imaging,” Appl. Opt., vol. 47, no. 10, pp. B44–B51, 2008. [10] J. M. Eichenholz and J. Dougherty, “Ultracompact fully integrated megapixel multispectral imager,” Integr. Opt., Devices, Mater., Technol., vol. 7218, pp. 721814-1–721814-10, Feb. 2009. [11] P. A. Stupar, R. L. Borwick, J. F. DeNatale, P. H. Kobrin, and W. J. Gunning, “Mems tunable fabry-perot filters with thick, two sided optical coatings,” in Proc. Int. Solid-State Sens., Actuators Microsyst. Conf., Jun. 2009, pp. 1357–1360. [12] D. Hays et al., “A hybrid mems–fiber optic tunable fabry–perot filter,” IEEE Microelectromech. Syst., J., vol. 19, no. 2, pp. 419–429, Apr. 2010. [13] A. Brückner, “Microoptical multi aperture imaging system,” Ph.D. dissertation, Dept. Comput. Sci., Fraunhofer Inst. Appl. Opt. Precision Eng., Jena, Germany, 2011. [14] S. Kim, E. Lee, M. H. Hayes, and J. Paik, “Multifocusing and depth estimation using a color shift model-based computational camera,” IEEE Trans. Image Process., vol. 21, no. 9, pp. 4152–4166, Sep. 2012. [15] Y. Bando, B. Chen, and T. Nishita, “Extracting depth and matte using a color-filtered aperture,” ACM Trans. Graph., vol. 27, no. 5, pp. 134:1–134:9, Dec. 2008. [16] S. Mendelson, A. Pajor, and N. Tomczak-Jaegermann, “Uniform uncertainty principle for bernoulli and subgaussian ensembles,” Construct. Approx., vol. 28, no. 3, pp. 277–289, 2008. [17] R. G. Antoninia, T. C. Hub, and A. Volodinc, “On the concentration phenomenon for φ-subgaussian random elements,” Statist. Probab. Lett., vol. 76, no. 5, pp. 465–469, Mar. 2006. [18] H. Arguello and G. R. Arce, “Rank minimization code aperture design for spectrally selective compressive imaging,” IEEE Trans. Image Process., vol. 22, no. 3, pp. 941–954, Mar. 2013. [19] E. Candes and T. Tao, “Decoding by linear programming,” IEEE Trans. Inf. Theory, vol. 51, no. 12, pp. 4203–4215, Dec. 2005. [20] H. Rauhut, “Compressive sensing and structured random matrices,” in Theoretical Foundations and Numerical Methods for Sparse Recovery, vol. 9, Berlin, Germany: Walter de Gruyter, 2010, pp. 1–92. [21] M. Davenport. (2011, Apr.). Proof of the RIP for Sub-Gaussian Matrices, Connexions [Online]. Available: http://cnx.org/content/m37186/1.4/ [22] U. Ayaz and H. Rauhut, “Nonuniform sparse recovery with subgaussian matrices,” Electron. Trans. Numerical Anal., vol. 32, no. 2, pp. 242–254, 2012. [23] M. Figueiredo, R. Nowak, and S. Wright, “Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems,” IEEE J. Sel. Topics Signal Process., vol. 1, no. 4, pp. 586–597, Dec. 2007.