# Undersampled Hyperspectral Image Reconstruction Based on Surfacelet Transform

Authors: Lei Liu; Jingwen Yan; Di Guo; Yunsong Liu; Xiaobo Qu
Journal: Journal of Sensors
 (2015)

---

## Abstract

Hyperspectral imaging is a crucial technique for military and environmental monitoring. However, limited equipment hardware resources severely affect the transmission and storage of a huge amount of data for hyperspectral images. This limitation has the potentials to be solved by compressive sensing (CS), which allows reconstructing images from undersampled measurements with low error. Sparsity and incoherence are two essential requirements for CS. In this paper, we introduce surfacelet, a directional multiresolution transform for 3D data, to sparsify the hyperspectral images. Besides, a Gram-Schmidt orthogonalization is used in CS random encoding matrix, two-dimensional and three-dimensional orthogonal CS random encoding matrixes and a patch-based CS encoding scheme are designed. The proposed surfacelet-based hyperspectral images reconstruction problem is solved by a fast iterative shrinkage-thresholding algorithm. Experiments demonstrate that reconstruction of spectral lines and spatial images is significantly improved using the proposed method than using conventional three-dimensional wavelets, and growing randomness of encoding matrix can further improve the quality of hyperspectral data. Patch-based CS encoding strategy can be used to deal with large data because data in different patches can be independently sampled.


---

## Body

## 1. Introduction

Typical hyperspectral imaging (HSI) is acquired on satellites and aerospace probes and then transmitted to grounds. The imaging spectrometer can provide tens to hundreds of narrow-band spectral information for each spatial location. The huge amounts of data but scarce equipment hardware resources on satellites and aerospace severely limit the transmission and storage of hyperspectral images [1, 2]. However, the ground receiving side holds very strong processing capability. If we transfer the system complexity from the satellites and aerospace probes to the ground, it is expected to potentially solve the limitation in traditional high spectral sampling.Traditionally, super-resolution reconstruction is used to improve the spatial resolution of HSI [3, 4] while compression can improve the transmission efficiency [5–8]. We acquire “all” data and then “throw away” most of it in transmission process. Can we just directly measure the part that we need to save the cost of storage? Compressed sensing (CS) has been proposed to solve these contradictions in remote sensing. The theory of CS shows that a sparse signal can be recovered from a relatively small number of linear measurements [9, 10]. The difference between the conventional method and CS based transmission method is shown in Figure 1.Figure 1
Conventional HSI transmission and CS HSI transmission. (a) Conventional sampling, compression, and decompression. (b) Compressed sensing and sparse reconstruction.
(a)
(b)CS has been widely used in medical imaging [11, 12] and wireless sensor network [13, 14] to reduce the sampled data in recent years which has been applied in remote sensing. Duarte et al. [15] proposed a single-pixel imaging method. Ma [2] applied CS in single frame imaging. Aravind et al. [16] used ten spectral bands to compare orthogonal matching pursuit with simultaneous orthogonal matching pursuit in reconstruction. Fowler [17] presented compressive-projection PCA to effectively shift the computational burden of PCA from the encoder to the decoder.According to the CS theory, the reconstruction error is bounded by the sparse approximation error [18]. For a signal x∈RN, let x~ be the solution to obey CS conditions, ΨT a forward transform, and (ΨTx)s the sparse approximation of x in transform domain by retaining only the s-largest entries of ΨTx; then, the reconstruction error is
(1)‖x~-x‖2≤C0ε+C1S‖ΨTx-(ΨTx)s‖1,
where C0, C1 are small positive constants and ε is noise power [19]. Equation () implies that a sparser representation will reduce the reconstruction error. An optimal sparsifying transform is always important for sparse image reconstruction to reduce the reconstruction error.In this study, we introduce the surfacelet transform (ST) to sparsely represent hyperspectral images by making use of the spatial and spectral information. The advantage of ST over wavelet on sparse HSI data reconstruction is demonstrated, where the data are reconstructed by using a fast iterative shrinkage-thresholding algorithm [20]. To further improve the reconstruction performance, a 3-dimensional (3D) random encoding is designed and the Gram-Schmidt orthogonalization is adopted. Finally, a patch-based CS encoding scheme is designed to deal with large size data.The remainder of this paper is organized as follows. First, ST-based compressive sensing HSI reconstruction method is introduced in Section2. The experimental results are presented in Section 3. Finally, conclusions are given in Section 4.


## 2. ST-Based Compressive Sensing HSI Reconstruction

The CS theory comes up with two principles: sparsity, which asserts that the reconstructed signal is sparse with a transformΨ, and incoherence, which requires that the encoding matrix Φ is incoherent with Ψ [10].


### 2.1. Sparse Representation of HSI Using Surfacelet Transform

Wavelet is commonly used as a typical sparse transform for HSI [21–24]. However, wavelet sometimes fails in sparsely representing HSI because it has only 7 directions [25–27]. Surfacelet transform (ST), proposed by Lu and Do [25], can efficiently capture the surface intrinsic geometrical structure within N-dimensional signals [25]. It offers 3×2L directional subbands with decomposition level L by combining the multiscale pyramid with the 3-dimensional directional filter banks (3D-DFB). Thus, surfacelet with more directions may help reducing the blocky artifacts caused by orthogonal wavelets [24, 25].ST is a multiscale version of the 3D-DFB. The input signalx[n] first goes through the 3D hourglass filter P0,00,0(ω), which is a three-channel undecimated filter bank. One branch of the three-channel structure of 3D-DFB is given by Figure 2. The output y[n] is then fed into a 2D filter bank, denoted by IRC12(l2), which operates on the (n1, n2) planes. The tree-structured filter bank IRC12(l2) produces 2l2 output subbands, denoted by zi[n] for 0≤i≤2l2. Each output is then fed into another 2D filter bank IRC13(l3) operating on the (n1, n3) planes. In the end, we get 2l2+l3 outputs, represented by zi,j[n] for 0≤i≤2l2 and 0≤j≤2l3 [25]. 3×2L spatial domain basis images of ST with 2 levels of decomposition are shown in Figure 3.Figure 2
One branch of the three-channel structure of 3D-DFB.Figure 3
3
×
2
L basis images of ST in 3D spatial domains with 2-level decomposition.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)ST is optimal for theC2 singularities, thus providing sparser representation of smooth curves and surface singularities than wavelet. Theoretically, for the C2 singularities of an image x, the best S-term approximation error ‖x~s-x‖2 using ST has error decay rate of O(S-2), while this error decay rate is O(S-1) for typical wavelets [28, 29]. Figure 4 shows that a sparser representation is achieved using ST. According to the CS theory [19], the reconstruction error bound is proportional to ‖ΨTx-(ΨTx)S‖1, where ΨT is the forward transform and S means the number of preserved coefficients in transform domain. Therefore, ST is expected to reduce the image reconstruction error of HSI.Figure 4
TheS-term approximation errors and reconstructed images when ‖ΨTx-(ΨTx)S‖1=4×103, (a) Lunar Lake HSI data, (b) the S-term approximation error, (c) the original 5th band, (d) sparsely represented image using wavelet (PSNR=33.12 dB), and (e) sparsely represented image using surfacelet (PSNR=38.46 dB).
(a)
(b)
(c)
(d)
(e)


### 2.2. Incoherence and Gaussian Random Encoding Matrix

Incoherence means a column of the CS encoding matrix must be trying the proliferation in the corresponding sparsity basisΨ. Since this paper focuses on investigating the spatial-spectral sparsity of HSI, a Gaussian random matrix is chosen as Φ because it is incoherent with the entire existing basis [9, 10]. If each spatial image is encoded separately, Φ can be represented by block matrix:
(2)Φ=[Φ10⋯00Φ2⋯0⋮⋱⋮00⋯ΦL],
where Φi(1≤i≤L) is an encoding matrix for the band i of HSI. When Φ1=Φ2=⋯=ΦL, that is, each 2D image is encoded with the same encoding matrix, then Φ is called 2D random encoding matrix; when Φ1≠Φ2≠⋯≠ΦL, that is, each 2D image is encoded differently, then Φ is called 3D random encoding matrix. Physically, 3D random encoding means that digital micromirror device (DMD) arrays are different for different bands. The performance of both encoding schemes will be discussed in Section 3.In order to improve the performance of the recovered signal, a Gram-Schmidt orthogonalization (GSO) is used in CS encoding matrix. Given one band random encoding matrixΦ-i(1≤i≤L) and the nth column denote as vn, the GSO is expressed by
(3)v-n=vn-∑m=1n-1〈vn,v-m〉vm,2≤n≤J,φn=v-n‖v-n‖22,1≤n≤J.The new matrixΦi=(φ1,φ2,…,φJ) is adopted as the CS encoding matrix of one spectral band. The randomness of the matrix is not affected by GSO. In the following, columns of each Φi are processed with GSO.


### 2.3. Reconstruction Algorithm

For HSIx=(x(1),x(2),…,x(L))∈RM×N×L, where x(i) represents the ith spectral band image, M×N represent the spatial dimensions, and L represents the spectral depth of HSI. Let xi=vec⁡(x(i))(1≤i≤L), where xi is a column vector of the ith band HSI with size MN×1. The data acquisition model for CS is given by
(4)y=[y1y2⋮yL]=[Φ10J×MN⋯0J×MN0J×MNΦ2⋯0J×MN⋮⋱⋮0J×MN0J×MN⋯ΦL][x1x2⋮xL],
where Φi∈RJ×MN(J<MN) is a random encoding matrix for the ith spectral band image and yi∈RJ×1 is the acquired undersampled data. The sampling ratio is defined as
(5)CSR=JMN,
where CSR<1, which means the HSI are undersampled.In this paper, ST is adopted to provide sparser representation of HSI data and is expected to reduce the reconstruction error. LetΨ represent inverse ST and let ΨTdenote forward ST; CS recovers x by solving [18, 19]
(6)α^=argmin⁡α⁡{λ‖α‖1+‖y-ΦΨα‖22},
where α=ΨTx, ‖·‖p(p=1,2) stands for lp-norm, and λ is the regularization parameter which decides the tradeoff between the sparsity and the data fidelity. Many researchers seek for a simple and fast algorithm to solve (), such as the conjugate gradient [30], Bregman iteration [31, 32], low rank reconstruction [33], and other methods. In this paper, we choose FISTA to solve () because of its simplicity and fast convergence, whose convergence rate is O(k-2), where k is an iteration counter [20].To solve (), the smallest Lipschitz constant [20] of the gradient ‖y-ΦΨα‖22 is G=2λmax⁡((ΦΨ)T(ΦΨ)). The Lipschitz constant ensures the convergence of algorithm [34]. Taking α0=0 and t1=1 as the initial values, the threshold is λ/G. For k≥1, solutions are found by iterating from () as follows:
(7)α^k=αk-1G(ΨTΦT(ΦΨαk-y)),αk=max⁡⁡{|α^k|-threshold,0}+sgn⁡(α^k),tk+1=1+1+4tk22,αk+1=αk+(tk-1tk+1)(αk-αk-1).
The maximum number of iterations of FISTA is set as 200 to achieve stable solutions. The final output is x^=Ψαk.


### 2.4. Patch-Based Compressed Sensing with Surfacelet Reconstruction (PCSST)

For a larger HSI datax, bigger encoding matrix may exceed the memory of computer or leads long computation time. In this case, a patch-based sampling operation can be used. Let R be a patch operator and let Rl(1≤l≤L) divide the lth band of x into n patches, and then the data acquisition model for CS is given by
(8)y=[y1y2⋮yL]=[Θ10⋯00Θ2⋯0⋮⋱⋮0⋯ΘL][R1x1R2x2⋮RLxL],
where yi(1≤i≤L) denotes measurements of patches in ith band of x and Θi satisfies Θi=(Θi1,Θi2,…,Θin)T, where Θij means random encoding on each block. In our scheme, the encoding on different patches of each band is different and these patches are nonoverlapped.PCSST allows implementing the proposed method on a larger dataset. It has the potential to allow for parallel computing as shown in Figure5 and also provides the possibility to reduce the complexity of sensor arrays since the sensing detectors independently sample the data in different patches.Figure 5
Parallel processing of compressive sampling on HIS.


## 2.1. Sparse Representation of HSI Using Surfacelet Transform

Wavelet is commonly used as a typical sparse transform for HSI [21–24]. However, wavelet sometimes fails in sparsely representing HSI because it has only 7 directions [25–27]. Surfacelet transform (ST), proposed by Lu and Do [25], can efficiently capture the surface intrinsic geometrical structure within N-dimensional signals [25]. It offers 3×2L directional subbands with decomposition level L by combining the multiscale pyramid with the 3-dimensional directional filter banks (3D-DFB). Thus, surfacelet with more directions may help reducing the blocky artifacts caused by orthogonal wavelets [24, 25].ST is a multiscale version of the 3D-DFB. The input signalx[n] first goes through the 3D hourglass filter P0,00,0(ω), which is a three-channel undecimated filter bank. One branch of the three-channel structure of 3D-DFB is given by Figure 2. The output y[n] is then fed into a 2D filter bank, denoted by IRC12(l2), which operates on the (n1, n2) planes. The tree-structured filter bank IRC12(l2) produces 2l2 output subbands, denoted by zi[n] for 0≤i≤2l2. Each output is then fed into another 2D filter bank IRC13(l3) operating on the (n1, n3) planes. In the end, we get 2l2+l3 outputs, represented by zi,j[n] for 0≤i≤2l2 and 0≤j≤2l3 [25]. 3×2L spatial domain basis images of ST with 2 levels of decomposition are shown in Figure 3.Figure 2
One branch of the three-channel structure of 3D-DFB.Figure 3
3
×
2
L basis images of ST in 3D spatial domains with 2-level decomposition.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)ST is optimal for theC2 singularities, thus providing sparser representation of smooth curves and surface singularities than wavelet. Theoretically, for the C2 singularities of an image x, the best S-term approximation error ‖x~s-x‖2 using ST has error decay rate of O(S-2), while this error decay rate is O(S-1) for typical wavelets [28, 29]. Figure 4 shows that a sparser representation is achieved using ST. According to the CS theory [19], the reconstruction error bound is proportional to ‖ΨTx-(ΨTx)S‖1, where ΨT is the forward transform and S means the number of preserved coefficients in transform domain. Therefore, ST is expected to reduce the image reconstruction error of HSI.Figure 4
TheS-term approximation errors and reconstructed images when ‖ΨTx-(ΨTx)S‖1=4×103, (a) Lunar Lake HSI data, (b) the S-term approximation error, (c) the original 5th band, (d) sparsely represented image using wavelet (PSNR=33.12 dB), and (e) sparsely represented image using surfacelet (PSNR=38.46 dB).
(a)
(b)
(c)
(d)
(e)


## 2.2. Incoherence and Gaussian Random Encoding Matrix

Incoherence means a column of the CS encoding matrix must be trying the proliferation in the corresponding sparsity basisΨ. Since this paper focuses on investigating the spatial-spectral sparsity of HSI, a Gaussian random matrix is chosen as Φ because it is incoherent with the entire existing basis [9, 10]. If each spatial image is encoded separately, Φ can be represented by block matrix:
(2)Φ=[Φ10⋯00Φ2⋯0⋮⋱⋮00⋯ΦL],
where Φi(1≤i≤L) is an encoding matrix for the band i of HSI. When Φ1=Φ2=⋯=ΦL, that is, each 2D image is encoded with the same encoding matrix, then Φ is called 2D random encoding matrix; when Φ1≠Φ2≠⋯≠ΦL, that is, each 2D image is encoded differently, then Φ is called 3D random encoding matrix. Physically, 3D random encoding means that digital micromirror device (DMD) arrays are different for different bands. The performance of both encoding schemes will be discussed in Section 3.In order to improve the performance of the recovered signal, a Gram-Schmidt orthogonalization (GSO) is used in CS encoding matrix. Given one band random encoding matrixΦ-i(1≤i≤L) and the nth column denote as vn, the GSO is expressed by
(3)v-n=vn-∑m=1n-1〈vn,v-m〉vm,2≤n≤J,φn=v-n‖v-n‖22,1≤n≤J.The new matrixΦi=(φ1,φ2,…,φJ) is adopted as the CS encoding matrix of one spectral band. The randomness of the matrix is not affected by GSO. In the following, columns of each Φi are processed with GSO.


## 2.3. Reconstruction Algorithm

For HSIx=(x(1),x(2),…,x(L))∈RM×N×L, where x(i) represents the ith spectral band image, M×N represent the spatial dimensions, and L represents the spectral depth of HSI. Let xi=vec⁡(x(i))(1≤i≤L), where xi is a column vector of the ith band HSI with size MN×1. The data acquisition model for CS is given by
(4)y=[y1y2⋮yL]=[Φ10J×MN⋯0J×MN0J×MNΦ2⋯0J×MN⋮⋱⋮0J×MN0J×MN⋯ΦL][x1x2⋮xL],
where Φi∈RJ×MN(J<MN) is a random encoding matrix for the ith spectral band image and yi∈RJ×1 is the acquired undersampled data. The sampling ratio is defined as
(5)CSR=JMN,
where CSR<1, which means the HSI are undersampled.In this paper, ST is adopted to provide sparser representation of HSI data and is expected to reduce the reconstruction error. LetΨ represent inverse ST and let ΨTdenote forward ST; CS recovers x by solving [18, 19]
(6)α^=argmin⁡α⁡{λ‖α‖1+‖y-ΦΨα‖22},
where α=ΨTx, ‖·‖p(p=1,2) stands for lp-norm, and λ is the regularization parameter which decides the tradeoff between the sparsity and the data fidelity. Many researchers seek for a simple and fast algorithm to solve (), such as the conjugate gradient [30], Bregman iteration [31, 32], low rank reconstruction [33], and other methods. In this paper, we choose FISTA to solve () because of its simplicity and fast convergence, whose convergence rate is O(k-2), where k is an iteration counter [20].To solve (), the smallest Lipschitz constant [20] of the gradient ‖y-ΦΨα‖22 is G=2λmax⁡((ΦΨ)T(ΦΨ)). The Lipschitz constant ensures the convergence of algorithm [34]. Taking α0=0 and t1=1 as the initial values, the threshold is λ/G. For k≥1, solutions are found by iterating from () as follows:
(7)α^k=αk-1G(ΨTΦT(ΦΨαk-y)),αk=max⁡⁡{|α^k|-threshold,0}+sgn⁡(α^k),tk+1=1+1+4tk22,αk+1=αk+(tk-1tk+1)(αk-αk-1).
The maximum number of iterations of FISTA is set as 200 to achieve stable solutions. The final output is x^=Ψαk.


## 2.4. Patch-Based Compressed Sensing with Surfacelet Reconstruction (PCSST)

For a larger HSI datax, bigger encoding matrix may exceed the memory of computer or leads long computation time. In this case, a patch-based sampling operation can be used. Let R be a patch operator and let Rl(1≤l≤L) divide the lth band of x into n patches, and then the data acquisition model for CS is given by
(8)y=[y1y2⋮yL]=[Θ10⋯00Θ2⋯0⋮⋱⋮0⋯ΘL][R1x1R2x2⋮RLxL],
where yi(1≤i≤L) denotes measurements of patches in ith band of x and Θi satisfies Θi=(Θi1,Θi2,…,Θin)T, where Θij means random encoding on each block. In our scheme, the encoding on different patches of each band is different and these patches are nonoverlapped.PCSST allows implementing the proposed method on a larger dataset. It has the potential to allow for parallel computing as shown in Figure5 and also provides the possibility to reduce the complexity of sensor arrays since the sensing detectors independently sample the data in different patches.Figure 5
Parallel processing of compressive sampling on HIS.


## 3. Experimental Results

Experiments are conducted in three aspects. First, better reconstructed signal using ST than using wavelet is demonstrated. Second, improving the randomness of encoding matrix for each spectral band is shown to improve the recovery quality and the reconstruction performance of two encoding schemes is compared. Third, a patch-based CS encoding scheme is designed to deal with large data.The HSI data is obtained from U.S. AVRIS website, including Moffet Field and Lunar Lake [35]. These data contain 224 spectral bands with spatial size 64×64, and every pixel is encoded with 16 bits. Linear interpolation is employed to fix the junk bands [36] and maintain the consistency of spectral lines, and all data are normalized.To evaluate the performance, the mean-square-error (MSE), peak signal-to-noise ratio (PSNR), and structural similarity (SSIM) [37] are adopted as criteria to measure the reconstruction error. Their definitions are
(9)MSE=1M·N∑j=1M∑i=1N(x^(i,j)-x(i,j))2,PSNR=10log10(1MSE),SSIM(x^,x)=(2μx^μx+C1)(2σx^x+C2)(μx^2+μx2+C1)(σx^2+σx2+C2),
where x denotes the original image, μx and σx are mean and standard deviation of x, x^ stands for the recovered image, μx^ and σx^ are mean and standard deviation of x^; C1, C2 are small constant, and C1=C2=0.08 is adopted in our work, which is used to avoid instability when either (μx^2+μx2) or (σx^2+σx2) is very close to zero [37]. Mean SSIM index is adopted in this paper. Simulations run on a dual core 2.5 GHz CPU laptop with 4 GB RAM. A fast ST implementation in C language is adopted.


### 3.1. Surfacelet and Wavelet with 2D Encoding

To simulate the CS data acquisition, the HSI are undersampled by random encoding matrixΦ. 50% sampled data means Φi∈R(642/2)×642for xi∈R642×1(1≤i≤224). In order to improve the quality of the recovered signal, Gram-Schmidt orthogonalization is used to obtain the CS encoding matrix. With the 2D encoding matrix (i.e., each band has the same encoding matrix), the ST-based reconstruction is compared with the 3D wavelets-based reconstruction. Daubechies filters “db4” with 2 levels of decomposition is used, and there are 7 directional subbands for 3D wavelet. The multiscale pyramid with 2 decomposition levels is chosen, and there are 12 directional subbands for ST. The complexity of ST is O(nlog10n) and wavelet is O(n).In order to well discuss the relationship betweenλ and PSNR, we compared the PSNR performance from λ=10 to λ=10-4 as shown in Figure 6. λ=10-2 is empirically chosen to give optimal PSNR in our work. As shown in Figure 7, ST significantly improves the PSNR of each spectral band more than 3D wavelet. Edges and curves are better reconstructed using ST than 3D wavelets as shown in Figure 8. At a given spatial location, the spectral line reconstructed using ST is much more consistent with the ground truth than wavelets as shown in Figure 9. Under different sampling ratios, ST achieves much better PSNRs than wavelet as shown in Figure 10, and more advantage of ST is seen at low ratios. Besides, ST shows better performance for the Lunar Lake dataset which has richer geometric textures than Moffet Field dataset. Running time (unit: second) based on ST and WT with different sampled rate is shown in Table 1.Table 1
Running time (unit: second) based on ST and WT with different sampled rate.
Sampling ratio
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
ST
458.9
629.7
833.8
1002.9
1224.2
1381.8
1544.9
1758.7
WT
989.8
1115.5
1359.7
1639.9
1871.6
2089.6
2324.1
2452.9
Note: the ST is implemented in C and WT is implemented in MATLAB.Figure 6
PSNR performance versus differentλ. Note: test is performed on Lunar Lake dataset at sampling ratio R=0.4.Figure 7
PSNR performance using wavelet and surfacelet when sampling ratio is 0.40. (a) Moffet Field dataset and (b) Lunar Lake dataset.
(a)
(b)Figure 8
The reconstructed 199th band of Moffet Field and 5th band of Lunar Lake using wavelets or surfacelet when sampling ratioR=0.40. (a) and (d) are original Moffet Field and Lunar Lake bands, (b) PSNR=24.31 dB, SSIM=0.9020 and (e) PSNR=32.91 dB, SSIM=0.9462 are recovered bands using wavelet, and (c) PSNR=30.75 dB, SSIM=0.9384 and (f) PSNR=35.88 dB, SSIM=0.9553 are reconstructed bands using surfacelet.
(a)
(b)
(c)
(d)
(e)
(f)Figure 9
Reconstructed spectral line when sampling ratioR=0.40. (a) Moffet Field dataset and (b) Lunar Lake dataset.
(a)
(b)Figure 10
PSNR and SSIM at different sampling ratios. (a), (c) PSNR and SSIM for Moffet Field dataset and (b), (d) PSNR and SSIM for Lunar Lake dataset.
(a)
(b)
(c)
(d)


### 3.2. 2D and 3D Encoding Matrix

The performance of 2D encoding and 3D encoding schemes is compared in this section. The recovered spatial images with two encoding schemes are shown in Figure11. Edges and curves are better reconstructed using 3D encoding than 2D encoding. The recovered spectral lines using 3D encoding are more consistent with ground truth than using 2D encoding as shown in Figure 12. The improvement is more obvious for Moffet Field dataset with textures in low intensities as shown in Figures 11(a) and 12(a). These results imply that increasing the encoding randomness among different bands will achieve better reconstruction.Figure 11
The reconstructed 199th band of Moffet Field and 5th band of Lunar Lake using 2D encoding and 3D encoding when sampling ratioR=0.3. (a) and (d) are original Moffet Field and Lunar Lake bands, (b) PSNR=26.84 dB, SSIM=0.9120 and (e) PSNR=33.65 dB, SSIM=0.9330 are recovered bands using 2D encoding matrix for each band and (c) PSNR=34.99 dB, SSIM=0.9600 and (f) PSNR=37.84 dB, SSIM=0.9701 are reconstructed bands using different encoding matrix for each band.
(a)
(b)
(c)
(d)
(e)
(f)Figure 12
Reconstructed spectral line when sampling ratioR=0.3. (a) Moffet Field dataset and (b) Lunar Lake dataset.
(a)
(b)


### 3.3. Reconstruction with PCSST

We tested the methods on Lunar Lake data with size256×256×224. A reconstructed spectral band is shown in Figure 13. The image structures recovered by surfacelet are much sharper than wavelet.Figure 13
Patch-based compressed sensing ST reconstruction. (a) The original 5th band of256×256 Lunar Lake dataset, (b) the patch-based reconstructed image using wavelet when undersampled rate R=0.4 and the PSNR=32.92 dB, SSIM=0.9551, and (c) the patch-based reconstructed image using ST when undersampled rate R=0.4 and the PSNR=34.99 dB, SSIM=0.9876.
(a)
(b)
(c)


## 3.1. Surfacelet and Wavelet with 2D Encoding

To simulate the CS data acquisition, the HSI are undersampled by random encoding matrixΦ. 50% sampled data means Φi∈R(642/2)×642for xi∈R642×1(1≤i≤224). In order to improve the quality of the recovered signal, Gram-Schmidt orthogonalization is used to obtain the CS encoding matrix. With the 2D encoding matrix (i.e., each band has the same encoding matrix), the ST-based reconstruction is compared with the 3D wavelets-based reconstruction. Daubechies filters “db4” with 2 levels of decomposition is used, and there are 7 directional subbands for 3D wavelet. The multiscale pyramid with 2 decomposition levels is chosen, and there are 12 directional subbands for ST. The complexity of ST is O(nlog10n) and wavelet is O(n).In order to well discuss the relationship betweenλ and PSNR, we compared the PSNR performance from λ=10 to λ=10-4 as shown in Figure 6. λ=10-2 is empirically chosen to give optimal PSNR in our work. As shown in Figure 7, ST significantly improves the PSNR of each spectral band more than 3D wavelet. Edges and curves are better reconstructed using ST than 3D wavelets as shown in Figure 8. At a given spatial location, the spectral line reconstructed using ST is much more consistent with the ground truth than wavelets as shown in Figure 9. Under different sampling ratios, ST achieves much better PSNRs than wavelet as shown in Figure 10, and more advantage of ST is seen at low ratios. Besides, ST shows better performance for the Lunar Lake dataset which has richer geometric textures than Moffet Field dataset. Running time (unit: second) based on ST and WT with different sampled rate is shown in Table 1.Table 1
Running time (unit: second) based on ST and WT with different sampled rate.
Sampling ratio
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
ST
458.9
629.7
833.8
1002.9
1224.2
1381.8
1544.9
1758.7
WT
989.8
1115.5
1359.7
1639.9
1871.6
2089.6
2324.1
2452.9
Note: the ST is implemented in C and WT is implemented in MATLAB.Figure 6
PSNR performance versus differentλ. Note: test is performed on Lunar Lake dataset at sampling ratio R=0.4.Figure 7
PSNR performance using wavelet and surfacelet when sampling ratio is 0.40. (a) Moffet Field dataset and (b) Lunar Lake dataset.
(a)
(b)Figure 8
The reconstructed 199th band of Moffet Field and 5th band of Lunar Lake using wavelets or surfacelet when sampling ratioR=0.40. (a) and (d) are original Moffet Field and Lunar Lake bands, (b) PSNR=24.31 dB, SSIM=0.9020 and (e) PSNR=32.91 dB, SSIM=0.9462 are recovered bands using wavelet, and (c) PSNR=30.75 dB, SSIM=0.9384 and (f) PSNR=35.88 dB, SSIM=0.9553 are reconstructed bands using surfacelet.
(a)
(b)
(c)
(d)
(e)
(f)Figure 9
Reconstructed spectral line when sampling ratioR=0.40. (a) Moffet Field dataset and (b) Lunar Lake dataset.
(a)
(b)Figure 10
PSNR and SSIM at different sampling ratios. (a), (c) PSNR and SSIM for Moffet Field dataset and (b), (d) PSNR and SSIM for Lunar Lake dataset.
(a)
(b)
(c)
(d)


## 3.2. 2D and 3D Encoding Matrix

The performance of 2D encoding and 3D encoding schemes is compared in this section. The recovered spatial images with two encoding schemes are shown in Figure11. Edges and curves are better reconstructed using 3D encoding than 2D encoding. The recovered spectral lines using 3D encoding are more consistent with ground truth than using 2D encoding as shown in Figure 12. The improvement is more obvious for Moffet Field dataset with textures in low intensities as shown in Figures 11(a) and 12(a). These results imply that increasing the encoding randomness among different bands will achieve better reconstruction.Figure 11
The reconstructed 199th band of Moffet Field and 5th band of Lunar Lake using 2D encoding and 3D encoding when sampling ratioR=0.3. (a) and (d) are original Moffet Field and Lunar Lake bands, (b) PSNR=26.84 dB, SSIM=0.9120 and (e) PSNR=33.65 dB, SSIM=0.9330 are recovered bands using 2D encoding matrix for each band and (c) PSNR=34.99 dB, SSIM=0.9600 and (f) PSNR=37.84 dB, SSIM=0.9701 are reconstructed bands using different encoding matrix for each band.
(a)
(b)
(c)
(d)
(e)
(f)Figure 12
Reconstructed spectral line when sampling ratioR=0.3. (a) Moffet Field dataset and (b) Lunar Lake dataset.
(a)
(b)


## 3.3. Reconstruction with PCSST

We tested the methods on Lunar Lake data with size256×256×224. A reconstructed spectral band is shown in Figure 13. The image structures recovered by surfacelet are much sharper than wavelet.Figure 13
Patch-based compressed sensing ST reconstruction. (a) The original 5th band of256×256 Lunar Lake dataset, (b) the patch-based reconstructed image using wavelet when undersampled rate R=0.4 and the PSNR=32.92 dB, SSIM=0.9551, and (c) the patch-based reconstructed image using ST when undersampled rate R=0.4 and the PSNR=34.99 dB, SSIM=0.9876.
(a)
(b)
(c)


## 4. Conclusions

Compressive sensing is a new sampling theorem. In this paper, the surfacelet transform is introduced into hyperspectral image reconstruction from compressive sampled data. The surfacelet transform is a directional multiresolution transform for 3D data, which is applied to sparsify the hyperspectral images for the first time. Simulations are conducted in three aspects. First, better reconstructed signal using surfacelet than using wavelet is demonstrated. Second, improving the randomness of encoding matrix for each spectral band is shown to improve the recovery quality. Third, a patch-based CS encoding scheme is designed to deal with large data. It provides the possibility to reduce the complexity of sensor arrays because sensing detectors in different patches independently sample data. Experiments demonstrate that reconstruction of spectral lines and spatial images is significantly improved using the proposed method than using conventional three-dimensional wavelets.In the future, our work includes the following two aspects.The first aspect is optimizing sparse representation of hyperspectral imaging, for example, combining adaptive sparse representation [38–40] and surfacelet transform. An adaptive dictionary may provide a sparser representation leading to a lower reconstruction error. Therefore, it is meaningful to try a low-complexity training method in the future.The second aspect is overlapping patches compressed sensing reconstruction. Overlapping patches can reduce the “block artifacts.” But overlapping also introduces more encoded data and higher computations on image reconstructions. Trading the compressive sampling rate with image qualities using overlapping and speeding up the reconstruction algorithm [41] will be discussed in the future.



---

*Source: 256391-2015-10-25.xml*