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

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 andMSSIM 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.


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][5][6][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.S ¸endur 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 algorithmhybrid 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][19][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.

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].

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. stands for the original image, whose size is  × .And then  is decomposed by binary wavelet transform with  scales, and we have where   is low-frequency component in the coarsest scale , while {  } =1,2,... is high-frequency component in different scale.Here,  = 1 is the finest scale.We now employ a sketch of the discrete curvelet transform algorithm.
(2) Let the size of initial block (the subband in the finest scale) be  min , and set  1 =  min .
( 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 ≈ length 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 where f  is the  biggest terms in the curvelet frame expansion to approximate .
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.

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  spectral bands and each band of HRSI has  ×  pixels.Let  denote a  ×  matrix of the  spectral observed vectors of size  ( =  × ).The  × 1 vector   ( = 1, 2, . . ., ) is the th column vector of the matrix , namely, lexicographically arranging th band image into a column vector   .In this paper, X is the vector predicted for the signal   of band .The  adjacent bands (not including itself) are utilized to perform MLR, where  is even; that is, where the  ×  matrix   consists of the adjacent column vectors of   (not including the th column vector).  has the form is the regression vector of size ×1.For  = 1, 2, . . ., , the least squares estimator of the regression vector   is given by β = (     ) −1      . (5)

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 (⋅) and  −1 (⋅) denote the 2-D curvelet transform and 2-D inverse curvelet transform.Let {  }  =1 be HRSI, where   is the th band image and  is the number of bands.  = (  ) stands for the coefficients of band   in the curvelet domain.A -level curvelet transform is performed on each   .Setting aside the coarsest coefficients    of   , the same scale and same directional coefficients from each high frequency coefficient    (   =   /   ) are stacked as a new 3-D datacube.The previous process is shown in Figure 1.In this way, several new 3-D datacubes {  }  =1 are obtained where  is the number of data.However, transforming {  }  =1 into   ( = 1, 2, . . .) is an inverse transform of the previous process.The spatial size of each datacube   varies with the scale and the direction, but each datacube has the same  spectral bands.The new three-dimensional datacubes are new representation of HRSI in the curvelet domain.The spectral correlation of   will be discussed.In this paper, each band   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,  = 25.
Let    (, ) denote the value of the position (, , ) of the datacube.Let  denote the set of spatial coordinates of   , and let || denote the element number of set .
The spectral correlation factor of   between band  and band  is defined as where 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   in the curvelet domain and the spectral correlation of the datacube Jasper Ridge.
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   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.

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 {  }  =1 for denoising.
The denoising process is as follows (Figure 3).

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 where   is the power of the pure signals (, , ), and   is the noise power in the noisy signals (, , ), while (, , ) stands for the position of the pixel in the HRSI; that is,   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.
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 × ] bands to get the best result through the lots of simulated experiments, the proposed algorithm utilizes [0.2 × ] bands to perform the MLR, where [] means getting the nearest even integer of .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.
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  denoised image is defined as and the MSSIM of band  denoised image is defined as where X  and    are the image contents at the th local window, and  is the number of local windows of the image.According to [30], the SSIM is defined as 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 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( ⋅ log( √ )) flops, where  is the number of pixels in the spatial domain and  is the number of bands.The complexity of the MLR is approximately the number of floating-point operations 4+2 2 + 3  flops where  and  are defined as before and  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.

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  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.

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.

Figure 2 :
Figure 2: The spectral correlation factor of {  }  =1 (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 3 :
Figure 3: Block diagram of the proposed method in this paper.

Figure 6 :
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.

Table 2 :
The computational time of AVIRIS data Cuprite, Jasper Ridge, Low Altitude, Lunar Lake, and Moffett Field (units: s).
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.Figure6shows band nos.20, 40, 60, 80, and 100 of the original image