Analysis and Denoising of Hyperspectral Remote Sensing Image in the Curvelet Domain

Abstract

A new denoising algorithm is proposed according to the characteristics of hyperspectral remote sensing image (HRSI) in the curvelet domain. Firstly, each band of HRSI is transformed into the curvelet domain, and the sets of subband images are obtained from different wavelength of HRSI. And then the detail subband images in the same scale and same direction from different wavelengths of HRSI are stacked to obtain new 3-D datacubes of the curvelet domain. Again, the characteristics analysis of these 3-D datacubes is performed. The analysis result shows that each new 3-D datacube has the strong spectral correlation. At last, due to the strong spectral correlation of new 3-D datacubes, the multiple linear regression is introduced to deal with these new 3-D datacubes in the curvelet domain. The simulated and the real data experiments are performed. The simulated data experimental results show that the proposed algorithm is superior to the compared algorithms in the references in terms of SNR. Furthermore, MSE and MSSIM in each band are utilized to show that the proposed algorithm is superior. The real data experimental results show that the proposed algorithm effectively removes the common spotty noise and the strip noise and simultaneously maintains more fine features during the denoising process.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Dong Xu 1 and Lei Sun 1 and Jianshu Luo 1 and Zhihui Liu 2

Recommended by Yue Wu

1, College of Science, National University of Defense Technology, Changsha, Hunan 410073, China
2, Department of Radiation Oncology, The 89th Hospital of PLA, Weifang, Shandong 261045, China

Received 7 May 2013; Accepted 1 July 2013

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Hyperspectral remote sensing image (HRSI) can be viewed as three-dimensional data consisting of one-dimensional spectral information and two-dimensional spatial information. With the fast development of hyperspectral remote sensing technology, HRSI can describe the characteristics of Earth objects more comprehensively and explicitly therefore, it is widely applied in many fields including agriculture, forestry, geological surveys, environmental monitoring, and military recon. Although over the last decades the development of imaging spectrometers is rapid, HRSI is still affected by many complex factors during the processing of acquisition and transmission, which will produce a mass of noises. The data that are contaminated with noise can cause a failure to extract valuable information and hamper further interpretation. In presence of noise in the image, extraction of all the useful information becomes difficult and noise can lead to artefacts and loss of spatial resolution [1]. The noise also affects the target detection, classification, and segmentation, so it is necessary to study the characteristic of HRSI for denoising [2].

Though the spatial correlation of HRSI is weaker than the nature image, the two-dimensional spatial information of HRSI is similar to the nature image [3]. When some traditional linear denoising methods are used to denoise the HRSI, they cannot get the satisfactory results for application such as Wiener filter. Compared with the nature image, the one-dimensional spectral information of HRSI is its particular characteristic. Minimum noise fraction (MNF), Savitzky-Golay filter, and wavelet denoising are the most popular among the existing denoising methods of hyperspectral imagery [4-7]. These methods are able to smooth the spectral while they have a negative impact on the sharp signal features.

In recent years, many denoising methods for HRSIs are constantly introduced with the development of hyperspectral remote sensing technology. Most of denoising methods combine correlation of spatial and spectral domain. Currently, wavelet denoising methods are widely used. The seminal work on signal denoising via wavelet thresholding or shrinkage proposed by Donoho and Johnstone [6, 7] shows that various wavelet thresholding schemes for denoising have near-optimal properties in the minimax sense. The observed noisy data in a local neighborhood is utilized to perform an approximate maximum a posteriori estimation of the variance for each coefficient [8]. Then, the predicted image coefficients are obtained by an approximate minimum mean squared error estimation procedure. Sendur and Selesnick developed a bivariate shrinkage function for image denoising [9, 10]. Their results show that the estimated wavelet coefficients depend on the parent coefficients. The smaller parent coefficients mean the greater shrinkage. In [11], Buades et al. proposed a new filter called the nonlocal means (NLM) filter, which takes account of the two most important attributes of a denoising algorithm: detail preservation and noise removal. Unfortunately, as most of filtering methods, NLM filter is also based on the additive Gaussian signal-independent (SI) noise assumption, so it cannot be directly applied for HRSI with signal-dependent (SD) noise [12]. In order to fully utilize the spatial information and the spectral information, Atkinson et al. proposed a denoising method that uses discrete Fourier transforms, 2-D discrete wavelet transforms, and soft thresholding of wavelet coefficients to denoise HRSI [13]. In [14], Othman and Qian developed a noise reduction algorithm--hybrid spatial-spectral noise-reduction algorithm (HSSNR) for hyperspectral datacube. The algorithm resorts to the spectral derivative domain, where the noise level is elevated, and benefits from the dissimilarity of the signal regularity in the spatial and spectral domains. Chen and Qian [15] combined wavelet denoising with dimensionality reduction for HRSI by using bivariate wavelet thresholding, wavelet packets, and principal component analysis (PCA). In [16], Chen and Qian developed a new denoising method using PCA and wavelet shrinkage to avoid removing the fine features of the datacubes during denoising process. PCA is used to decorrelate the fine features of the datacube from noise for detail preservation. In this paper the method in [16] is named as PCABS.

Since wavelet has good time-frequency-localization property and multiresolution analysis property, it is widely and successfully applied in several fields [17]. However, wavelet is not perfect. Wavelet is mainly applied to representation of isotropic singularity object, while for anisotropic singularity, such as boundary and linear features of an image, it is not a good representation tool. In other words wavelet is good representation for point singularity, but for the curves it is relatively weak [18-20]. The complex wavelet [18], curvelet transform [19], and contourlet transform [20] are proposed to overcome the drawbacks of the wavelet transform. This paper focuses on the curvelet transform applied for HRSI.

Compared with normal three-dimensional data cube of fixed variance of additive noise, the noise level of HRSI may vary dramatically from band to band. The noise standard deviation in each band of HRSI is not constant; in particular, there exist some bands at which the atmosphere absorbs so much light that the signal received from the surface is unreliable [21]. And SD noise may no longer be neglected with respect to the SI noise in the new-generation hyperspectral sensors [22, 23]. In this paper, the proposed denoising method is available for the both SD noise and SI noise. The representation of HRSI in another domain is proposed for detail preservation. By analyzing spectral correlation of HRSI in the new representation domain, we know the fact that the spectral correlation of HRSI in the new representation domain is strong. So the multiple linear regression (MLR) method is proposed to predict the representation of pure signal in the curvelet domain, in order to remove the noise of HRSI.

In the rest of this paper, the mathematical tools curvelet transform and MLR model are introduced in Section 2. Section 3 gives the analysis of HRSI in the curvelet domain. In Section 4 the denoising process is proposed in the curvelet domain. The simulated data experiment and the real data experiment are performed, and the experimental results are presented in Section 5. Finally, Section 6 draws the conclusion.

2. Mathematical Tools

In the section, the mathematical tools curvelet transform and MLR model are introduced. The curvelet transform, pioneered by Candès and Donoho, is shown to be optimal in a certain sense for functions in the domain with curved singularities [19]. Due to the strong spectral correlation of HRSI, MLR model has been widely applied to the HRSI [24, 25].

2.1. The Curvelet Transform

The curvelet transform is a new multiresolution analysis framework and widely applied in various image processing problems. The curvelet decomposition can be equivalently stated in the following four steps: (1) subband decomposition, (2) smooth partitioning, (3) renormalization, and (4) ridgelet transform. In short, the curvelet obtained by bandpass filtering of multiscale ridgelets with passband is rigidly linked to the scale of spatial localization. The discrete curvelet transform [26] is given as follows.

The digital image is used as an example to introduce the discrete curvelet transform. I stands for the original image, whose size is N × N . And then I is decomposed by binary wavelet transform with J scales, and we have [figure omitted; refer to PDF] where _{C J} is low-frequency component in the coarsest scale J , while _{{ D j } j = 1,2 , ... J} is high-frequency component in different scale. Here, j = 1 is the finest scale. We now employ a sketch of the discrete curvelet transform algorithm.

(1) Perform the binary wavelet transform with J scales on the original image I , and obtain subband sequence _{C J} and _{{ D j } j = 1,2 , ... J} .

(2) Let the size of initial block (the subband in the finest scale) be _{B min} , and set _{B 1} = _{B min} .

(3) For j = 1,2 , ... , J

(a) partition the subband _{D j} with a block size _{B j} ;

(b) apply the digital ridgelet transform to each block;

(c) If j module 2 = 1 _{B j + 1} = 2 _{B j} , and else _{B j + 1} = _{B j} .

Since each step of the previous decomposition process is invertible, the inverse curvelet transform is an invertible process.

The curvelet transform overcomes the major drawback that wavelets cannot really represent two-dimensional objects with edges sparsely and captures more directional information besides the horizontal, vertical, and diagonal directions. The system approximately obeys the scale relationship width [approximate] lengt ^{h 2} . Therefore, it is optimal in a certain sense for functions in the domain with curved singularities.

It is known that the curvelet transform competes surprisingly well with the ideal adaptive rate. The approximation error is obtained as [figure omitted; refer to PDF] where ^{f ~ m C} is the m biggest terms in the curvelet frame expansion to approximate f .

Curvelet is optimal in the sense that no other representation can yield a smaller asymptotic error with the same number of terms. Because of its surprising properties for image processing, a fast and accurate discrete curvelet transform operating on digital data is necessary. Candès et al. [27, 28] presented two 2-D discrete curvelet transforms for the second generation curvelet, which is curvelet via transforms and curvelet via wrapping of specially selected Fourier samples. Compared with the first generation curvelet, they are conceptually simpler, faster, and far less redundant. The curvelet transform used in this work is based on unequally-spaced fast Fourier transform.

2.2. Multiple Linear Regression (MLR)

Due to the strong spectral correlation of new 3-D datacubes in the curvelet domain, the MLR model is introduced to predict the representation of pure HRSI in the curvelet domain. It is assumed that the HRSI has B spectral bands and each band of HRSI has M × N pixels. Let X denote a P × B matrix of the B spectral observed vectors of size P ( P = M × N ) . The P × 1 vector _{X k} ( k = 1,2 , ... , B ) is the k th column vector of the matrix X , namely, lexicographically arranging k th band image into a column vector _{X k} . In this paper, _{X ^ k} is the vector predicted for the signal _{X k} of band k . The L adjacent bands (not including itself) are utilized to perform MLR, where L is even; that is, [figure omitted; refer to PDF] where the P × L matrix _{X λ k} consists of the adjacent column vectors of _{X k} (not including the k th column vector). _{X λ k} has the form [figure omitted; refer to PDF]

_{β k} is the regression vector of size L × 1 . For k = 1,2 , ... , B , the least squares estimator of the regression vector _{β k} is given by [figure omitted; refer to PDF]

3. Analysis of HRSI in the Curvelet Domain

HRSI is a datacube, having two spatial dimensions and a third spectral dimension. Fixing the wavelength band yields a 2-D image of the scene at a particular wavelength. So it may also be visualized as a stack of 2-D band images, each corresponding to a certain wavelength [3]. Since the bands are so closely spaced in wavelength, images in adjacent bands are highly correlated. In the rest of the section, HRSI is transformed into the curvelet domain band by band and the spectral correlation of the HRSI in the curvelet domain will be mainly discussed to denoise by performing the MLR.

Let C ( · ) and ^{C - 1} ( · ) denote the 2-D curvelet transform and 2-D inverse curvelet transform. Let { _{X i} ^{} i = 1 B} be HRSI, where _{X i} is the i th band image and B is the number of bands. _{C i} = C ( _{X i} ) stands for the coefficients of band _{X i} in the curvelet domain. A J -level curvelet transform is performed on each _{X i} . Setting aside the coarsest coefficients ^{C i L} of _{C i} , the same scale and same directional coefficients from each high frequency coefficient ^{C i H} ( ^{C i H} = _{C i} / ^{C i L} ) are stacked as a new 3-D datacube. The previous process is shown in Figure 1. In this way, several new 3-D datacubes { ^{F t } t = 1 S} are obtained where S is the number of data. However, transforming { ^{F t } t = 1 S} into _{C i} ( i = 1,2 , ... B ) is an inverse transform of the previous process. The spatial size of each datacube ^{F t} varies with the scale and the direction, but each datacube has the same B spectral bands. The new three-dimensional datacubes are new representation of HRSI in the curvelet domain. The spectral correlation of ^{F t} will be discussed. In this paper, each band _{X i} is decomposed into 5 levels. The first scale is coarsest scale with only one direction. From the second scale to the fourth scale, there are eight directions. At the finest scale, one direction is obtained. Setting aside the first scale, the total number of the other directions from different scale is 25; that is, S = 25 .

Figure 1: The process of obtaining new 3-D datacubes { ^{F t } t = 1 S} .

[figure omitted; refer to PDF]

Let ^{F l t} ( m , n ) denote the value of the position ( m , n , l ) of the datacube. Let H denote the set of spatial coordinates of ^{F t} , and let | H | denote the element number of set H .

The spectral correlation factor of ^{F t} between band k and band l is defined as [figure omitted; refer to PDF] where _{u k} = ( 1 / | H | ) _{∑ ( m , n ) ∈ H} ^{F k t} ( m , n ) and _{u l} = ( 1 / | H | ) _{∑ ( m , n ) ∈ H} ^{F l t} ( m , n ) . In this paper, let k = l + 1 ; thus _{ρ t} ( k , k + 1 ) stands for the correlation factor of neighboring bands for datacube ^{F t} . In order to analyze the spectral correlation of HRSI in the curvelet domain, AVIRIS hyperspectral data Jasper Ridge whose size is 256 × 256 × 224 (width × height × band) is used as an experimental example. Figure 2 shows the spectral correlation of each datacube ^{F t} in the curvelet domain and the spectral correlation of the datacube Jasper Ridge.

Figure 2: The spectral correlation factor of { ^{F t } t = 1 S} (a) eight directions of the second scale, (b) eight directions of the third scale, (c) eight directions of the fourth scale, (d) the finest scale, and (e) the spectral correlation factor of the Jasper Ridge.

[figure omitted; refer to PDF]

From Figure 2, it is clear that the absolute value of linear correlation factor between neighboring bands is close to 1 from the second scale to the fourth scale except for the bands 1-5, 105-115, and 150-170, which is due to water absorption and low signal-to-noise ratio (SNR) for the HRSI. These bands are often viewed as junk band. From the second scale to the finest scale, the correlation factor of junk bands in the curvelet domain becomes much lower, but the other bands are still keeping close to 1. The result in the curvelet domain is consistent with the original HRSI. It also indicates that the curvelet coefficients of hyperspectral data between neighboring bands have significant linear correlation from the same scale and the same direction. The curvelet transform keeps the strong spectral correlation of HRSI, and specially at some directions from the second scale to the forth scale the spectral correlation in the curvelet domain becomes stronger. The reason for this phenomenon is the strong spectral correlation of HRSI, which means that the pixels in the same spatial location of each band image are similar. Images at different wavelengths capture the same scene, which has the same physical structure; thus the profiles of different wavelength images are similar [29]. After the curvelet transform, the coefficients of datacube ^{F t} at close wavelengths maintain this similarity as well and have significant linear correlation. Since the most of subtle noise of HRSI is maintained in the finest scale, the correlation factor of the finest scale is lower than the other scales. But in the finest scale, the most of the correlation factors are still greater than 0.85. Due to the strong spectral correlation of HRSI in the curvelet domain, a new denoising algorithm based on MLR is proposed in the next section.

4. Proposed Denoising Algorithm

In this section, we summarize our denoising algorithm. According to the correlation factor of HRSI in the curvelet domain in Section 3, the MLR model is performed on the datacubes { ^{F t } t = 1 S} for denoising.

The denoising process is as follows (Figure 3).

(1) Input the noisy datacube X .

(2) A J -level curvelet transform (CT) is performed on _{X i} , and _{C i} = C ( _{X i} ) ( i = 1,2 , ... B ) is obtained.

(3) Set aside the coarsest scale coefficients ^{C i L} of _{C i} , and then the same scale and same directional coefficients from each high frequency coefficients ^{C i H} ( ^{C i H} = _{C i} / ^{C i L} ) are stacked as a 3-D datacube. The 3-D datacubes { ^{F t } t = 1 S} are obtained.

(4) MLR is performed on ^{F t} , t = 1,2 , ... S , and ^{F ^ t} ( t = 1,2 , ... S ) is obtained.

(5) For i = 1,2 , ... B , combine the new high frequency coefficients ^{C ^ i H} extracted from each ^{F ^ i t} ( t = 1,2 , ... S ) with the corresponding coarsest scale coefficients ^{C i L} to obtain the prediction coefficients _{C ^ i} of the _{C i} .

(6) The inverse curvelet transform (ICT) is performed to obtain the denoised image _{X ^ i} ; that is, _{X ^ i} = ^{C - 1} ( _{C ^ i} ) .

(7) Output the denoised datacube X ^ .

Figure 3: Block diagram of the proposed method in this paper.

[figure omitted; refer to PDF]

5. Experiment

Signal-to-noise ratio (SNR) is a key parameter on measuring the HRSI quality. So in this paper, we utilize SNR to evaluate the proposed algorithm. Here the SNR is defined as [figure omitted; refer to PDF] where _{P X} is the power of the pure signals P X ( i , j , k ) , and _{P N} is the noise power in the noisy signals X ( i , j , k ) , while ( i , j , k ) stands for the position of the pixel in the HRSI; that is, [figure omitted; refer to PDF]

5.1. Simulated Data Experiment

The simulated experiment of the noise reduction is carried out on AVIRIS images, Cuprite, Jasper Ridge, Low Altitude, Lunar Lake, and Moffett Field provided by JPL, NASA. The size of datacube we extracted from the Cuprite, Jasper Ridge, Low Altitude, Lunar Lake, and Moffett Field for testing is 256 × 256 × 224 (width × height × band). Figure 4 shows the band no.40 of these HRSIs. For the AVIRIS images a 28 m × 28 m ground sample distance (GSD) datacube is derived by spatially averaging the 4 m × 4 m GSD datacube elevating the nominal SNR to 7000 : 1. Having such high SNR, this datacube is viewed as a pure datacube [14], which is used as a reference to measure the SNR before and after denoising. The image is corrupted by Gaussian white noise. It is different from the simple stationary additive noise model that is simulated by adding noise with a fixed standard deviation to the datacube; the noise variance is proportional to the average signal amplitude of each band, but the noise in each band is still additive noise. The SNR of the simulated noisy data is 600 : 1, which is chosen by comprehensive requirement of the users and the machine design parameters.

Band no.80 of AVIRIS images.

(a) Cuprite

[figure omitted; refer to PDF]

(b) Jasper Ridge

[figure omitted; refer to PDF]

(c) Low Altitude

[figure omitted; refer to PDF]

(d) Lunar Lake

[figure omitted; refer to PDF]

(e) Moffett Field

[figure omitted; refer to PDF]

In order to indicate that the nature image denoising methods cannot be immediately used for HRSI noise reduction, the 2-D complex wavelet with bivariate shrinkage (CWBS) [10] and the curvelet denoising (CD) [26] are used for comparison. From Table 1, it is shown that the denoised results in terms of SNR are obtained by the two nature image denoising methods; even the lower SNR is obtained. In order to illustrate the superiority of the proposed algorithm in our paper, the HSSNR [14] and the PCABS [16] are also performed to compare with the proposed method in this paper. The parameters chosen for the PCABS are consistent with [16]. In this paper, through utilizing [ 0.2 × B ] bands to get the best result through the lots of simulated experiments, the proposed algorithm utilizes [ 0.2 × B ] bands to perform the MLR, where [ x ] means getting the nearest even integer of x . Table 1 shows the SNR of the hyperspectral datacubes Cuprite, Jasper Ridge, Low Altitude, Lunar Lake, and Moffett Field after denoising by the HSSNR, the PCABS, and the proposed algorithm. The results indicate that the proposed method is best method for denoising HRSI in terms of SNR.

Table 1: SNR of AVIRIS data Cuprite, Jasper Ridge, Low Altitude, Lunar Lake, and Moffett Field.

	SNR	CWBS	CD	HSSNR	PCABS	The proposed
Cuprite	600	1368.3	907.6	2712.6	8192.8	10952.0
Jasper Ridge	720.6	387.8	2365.9	6135.2	7503.0
Lunar Lake	2869.5	1844.2	2791.4	9787.5	12516.0
Low Altitude	530.7	220.8	2000.4	4351.5	5140.3
Moffett Field	530.1	253.7	2267.8	5757.7	6427.1

In order to deeply analyze the proposed algorithm, the mean square of errors (MSE) in each band and the mean structural similarity (MSSIM) [30] between before and after denoising are utilized in each band. The lower MSE means that the denoised image is more similar to the original pure image, while MSSIM is close to 1 means that the denoised image is more similar to the original pure image. The MSE of band k denoised image is defined as [figure omitted; refer to PDF] and the MSSIM of band k denoised image is defined as [figure omitted; refer to PDF] where ^{X ^ t k} and P ^{X t k} are the image contents at the q th local window, and Q is the number of local windows of the image. According to [30], the SSIM is defined as [figure omitted; refer to PDF] where _{C 1} = ( _{K 1} V ^{) 2} and _{C 2} = ( _{K 2} V ^{) 2} ( _{K 1} ...a; 1 and _{K 2} ...a; 1 are two small constants, _{K 1} = 0.01 and _{K 2} = 0.03 ) and V is the dynamic range of the pixel values. In this paper, we use 11 × 11 circular-symmetric Gaussian weighting function w = { _{w i} |" i = 1,2 , ... , N } , with standard deviation of 1.5 samples, normalized to unit sum ( ^{∑ i = 1 N} _{w i} = 1 ). The estimates of local statistics _{μ x} , _{σ x} and _{σ x y} are then modified accordingly as [figure omitted; refer to PDF] where _{μ y} and _{σ y} are consistent with the form _{μ x} and _{σ x} .

From Table 1, it is shown that the PCABS method and the proposed method obtain higher SNR. As a result, Figure 5 only shows MSE and MSSIM in each band of the AVIRIS datacubes denoised by the proposed algorithm and the PCABS. The MSE of each band obtained by the proposed algorithm is lower than the PCABS in nearly all bands, and the MSSIM of each band obtained by the proposed algorithm is closer to 1 than the PCABS, which indicates that the proposed algorithm is superior to PCABS.

The MSE and the MSSIM in each band obtained by PCABS and the proposed method.

(a) Cuprite

[figure omitted; refer to PDF]

(b) Jasper Ridge

[figure omitted; refer to PDF]

(c) Lunar Lake

[figure omitted; refer to PDF]

(d) Low Altitude

[figure omitted; refer to PDF]

(e) Moffett Field

[figure omitted; refer to PDF]

The computational complexity, that is, the number of floating-point operations (flops), of the proposed method can be analyzed as follows. The complexity of the curvelet transform is in order of O ( P B · log ( P ) ) flops, where P is the number of pixels in the spatial domain and B is the number of bands. The complexity of the MLR is approximately the number of floating-point operations 4 P L B + 2 P ^{L 2} B + ^{L 3} B flops where P and B are defined as before and L is the number of bands used to perform MLR. Therefore, the computational complexity of the proposed method is mainly contributed by the MLR. The computational complexity of the proposed method is greater than that of the HSSNR and PCABS. The computational times of the HSSNR, the PCABS, and the proposed method are given in Table 2. The previous denoising approaches are implemented using the Matlab programming language and run on a PC with a Pentium 2.70 GHz Dual-Core CPU and a 1024 MB RAM.

Table 2: The computational time of AVIRIS data Cuprite, Jasper Ridge, Low Altitude, Lunar Lake, and Moffett Field (units: s).

	HSSNR	PCABS	The proposed
Cuprite	145.12	185.25	810.09
Jasper Ridge	145.12	188.29	847.60
Lunar Lake	148.58	184.10	809.97
Low Altitude	146.98	201.05	860.91
Moffett Field	146.65	185.30	858.04

5.2. Real Data Experiment

In this paper the OMIS (operational modular imaging spectrometer) data that is developed by the Shanghai Institute of Technical Physics of the Chinese Academy of Sciences is used for real data experiment to verify the correctness and performance of algorithm. It has 128 spectral bands ranging from visible to thermal infrared wavelength. The size of datacube we extracted from OMIS data for testing is 256 × 256 × 128 (width × height × band). We perform denoising for the original OMIS data. Figure 6 shows band nos. 20, 40, 60, 80, and 100 of the original image and the denoised image obtained by PCABS and the proposed method, as well as the difference between the original image and the denoised image obtained by PCABS and the proposed method. In Figure 6, the denoised images (b) and (c) show that the proposed method is feasible; especially the band no. 80 can be seen obviously. The image (d) is the difference between the original image and the denoised image obtained by PCABS, and the image (e) is the difference between the original image and the denoised image obtained by the proposed algorithm. The difference shows that the two algorithms can effectively remove the common spotty noise and the strip noise from the noisy OMIS. But the features such as the edge of the image can be observed in the difference. The reason for this phenomenon is that the fine features are removed as noise during the denoising process. So the less fine features in the difference mean the better denoised results. From Figures 6(d) and 6(e), the image (d) contains more fine features than the image (e), which is easy to distinguish by human eyes. The result shows that the PCABS method removes more fine features during the denoising process. Therefore the proposed algorithm is superior to the PCABS in terms of removing noise and simultaneously maintaining fine features during the denoising process.

Figure 6: Band nos. 20, 40, 60, 80, and 100 of (a) the original image and (b) the denoised image obtained by the PCABS, (c) the proposed algorithm, and (d) the difference between the original image and the denoised image obtained by PCABS, as well as (e) the difference between the original image and the denoised image obtained by the proposed algorithm.

[figure omitted; refer to PDF]

6. Conclusions

In this paper, the spectral correlation of HRSI in the curvelet domain is discussed. By the analysis, in the curvelet domain, the strong spectral correlation of the hyperspectral remote sensing image is kept; even in some directions and scales it becomes stronger. So a new denoising algorithm is proposed; the MLR is performed in the curvelet domain to denoise the HRSI.

Simulated experimental result shows that the proposed method improves the quality of HRSI significantly in terms of SNR, MSE of each band, and MSSIM of each band. It is also seen that the denoised results obtained by the two algorithms are not content in the bands 1-5. This is a problem that will be studied in the future. For the real OMIS data, the results show that the proposed method is valid. The proposed method obtains better results in terms of detail preservation and noise removal during the denoising process.

Acknowledgment

This work was supported by the National Natural Science Foundation of China under Project 61101183 and Project 41201363.

References

[1] D. Letexier, S. Bourennane, "Noise removal from hyperspectral images by multidimensional filtering," IEEE Transactions on Geoscience and Remote Sensing , vol. 46, no. 7, pp. 2061-2069, 2008.

[2] L. Ebadi, H. M. Shafri, S. Mansor, R. Ashurov, "A review of applying second-generation wavelets for noise removal from remote sensing data," Environmental Earth Sciences , pp. 1-12, 2013.

[3] A. C. Zelinski, V. K. Goyal, "Denoising hyperspectral imagery and recovering junk bands using wavelets and sparse approximation," in Proceedings of IEEE International Geoscience and Remote Sensing Symposium (IGARSS 06), pp. 387-390, August 2006.

[4] A. A. Green, M. Berman, P. Switzer, M. D. Craig, "Transformation for ordering multispectral data in terms of image quality with implications for noise removal," IEEE Transactions on Geoscience and Remote Sensing , vol. 26, no. 1, pp. 65-74, 1988.

[5] J. Chen, P. Jönsson, M. Tamura, Z. Gu, B. Matsushita, L. Eklundh, "A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter," Remote Sensing of Environment , vol. 91, no. 3-4, pp. 332-344, 2004.

[6] D. L. Donoho, J. M. Johnstone, "Ideal spatial adaptation by wavelet shrinkage," Biometrika , vol. 81, no. 3, pp. 425-455, 1994.

[7] D. L. Donoho, "De-noising by soft-thresholding," IEEE Transactions on Information Theory , vol. 41, no. 3, pp. 613-627, 1995.

[8] M. K. Mihçak, I. Kozintsev, K. Ramchandran, P. Moulin, "Low-complexity image denoising based on statistical modeling of wavelet coefficients," IEEE Signal Processing Letters , vol. 6, no. 12, pp. 300-303, 1999.

[9] L. Sendur, I. W. Selesnick, "Bivariate shrinkage functions for wavelet-based denoising exploiting interscale dependency," IEEE Transactions on Signal Processing , vol. 50, no. 11, pp. 2744-2756, 2002.

[10] L. Sendur, I. W. Selesnick, "Bivariate shrinkage with local variance estimation," IEEE Signal Processing Letters , vol. 9, no. 12, pp. 438-441, 2002.

[11] A. Buades, B. Coll, J. M. Morel, "A review of image denoising algorithms, with a new one," Multiscale Modeling and Simulation , vol. 4, no. 2, pp. 490-530, 2005.

[12] A. Buades, B. Coll, J.-M. Morel, "A non-local algorithm for image denoising," in Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR '05), vol. 2, pp. 60-65, June 2005.

[13] I. Atkinson, F. Kamalabadi, D. L. Jones, "Wavelet-based hyperspectral image estimation," in IEEE International Geoscience and Remote Sensing Symposium (IGARSS '03), pp. 743-745, July 2003.

[14] H. Othman, S.-E. Qian, "Noise reduction of hyperspectral imagery using hybrid spatial-spectral derivative-domain wavelet shrinkage," IEEE Transactions on Geoscience and Remote Sensing , vol. 44, no. 2, pp. 397-408, 2006.

[15] G. Chen, S.-E. Qian, "Denoising and dimensionality reduction of hyperspectral imagery using wavelet packets, neighbour shrinking and principal component analysis," International Journal of Remote Sensing , vol. 30, no. 18, pp. 4889-4895, 2009.

[16] G. Chen, S.-E. Qian, "Denoising of hyperspectral imagery using principal component analysis and wavelet shrinkage," IEEE Transactions on Geoscience and Remote Sensing , vol. 49, no. 3, pp. 973-980, 2011.

[17] T. H. Yi, H. N. Li, H. N. Li, X. Y. Zhao, "Noise smoothing for structural vibration test signals using an improved wavelet thresholding technique," Sensors , vol. 12, pp. 11205-11220, 2012.

[18] I. W. Selesnick, R. G. Baraniuk, N. G. Kingsbury, "The dual-tree complex wavelet transform," IEEE Signal Processing Magazine , vol. 22, no. 6, pp. 123-151, 2005.

[19] E. J. Candes, D. L. Donoho, "Curvelets--a surprisingly effective nonadaptive representation for objects with edges," in Curve and Surface Fitting, pp. 105-120, Saint-Malo, France, 2000.

[20] M. N. Do, M. Vetterli, "The contourlet transform: an efficient directional multiresolution image representation," IEEE Transactions on Image Processing , vol. 14, no. 12, pp. 2091-2106, 2005.

[21] J. A. Richards, X. Jia Remote Sensing Digital Image Analysis , Springer, Berlin, Germany, 1999., 3rd.

[22] M. L. Uss, B. Vozel, V. V. Lukin, K. Chehdi, "Local signal-dependent noise variance estimation from hyperspectral textural images," IEEE Journal on Selected Topics in Signal Processing , vol. 5, no. 3, pp. 469-486, 2011.

[23] N. Acito, M. Diani, G. Corsini, "Signal-dependent noise modeling and model parameter estimation in hyperspectral images," IEEE Transactions on Geoscience and Remote Sensing , vol. 49, no. 8, pp. 2957-2971, 2011.

[24] J. M. Bioucas-Dias, J. M. P. Nascimento, "Hyperspectral subspace identification," IEEE Transactions on Geoscience and Remote Sensing , vol. 46, no. 8, pp. 2435-2445, 2008.

[25] N. Acito, M. Diani, G. Corsini, "Residual striping reduction in hyperspectral images," in Proceedings of the 17th International Conference on Digital Signal Processing (DSP '11), pp. 1-7, July 2011.

[26] J.-L. Starck, E. J. Candès, D. L. Donoho, "The curvelet transform for image denoising," IEEE Transactions on Image Processing , vol. 11, no. 6, pp. 670-684, 2002.

[27] E. J. Candès, D. L. Donoho, "New tight frames of curvelets and optimal representations of objects with piecewise C² singularities," Communications on Pure and Applied Mathematics , vol. 57, no. 2, pp. 0219-0266, 2004.

[28] E. Candès, L. Demanet, D. Donoho, L. Ying, "Fast discrete curvelet transforms," Multiscale Modeling and Simulation , vol. 5, no. 3, pp. 861-899, 2006.

[29] L. Sun, J.-S. Luo, "Junk band recovery for hyperspectral image based on curvelet transform," Journal of Central South University of Technology , vol. 18, no. 3, pp. 816-822, 2011.

[30] Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli, "Image quality assessment: from error visibility to structural similarity," IEEE Transactions on Image Processing , vol. 13, no. 4, pp. 600-612, 2004.