Matlab Code for: "Non-local Meets Global: An Integrated Paradigm for Hyperspectral Denoising. Arvix. Dec 2018"
Non-local low-rank tensor approximation has been developed as a state-of-the-art method for hyperspectral image (HSI) denoising. Unfortunately, with more spectral bands for HSI, while the running time of these methods significantly increases, their denoising performance benefits little. In this paper, we claim that the HSI underlines a global spectral low-rank subspace, and the spectral subspaces of each full band patch groups should underlie this global low-rank subspace. This motivates us to propose a unified spatial-spectral paradigm for HSI denoising. As the new model is hard to optimize, we further propose an efficient algorithm for optimization, which is motivated by alternating minimization. This is done by first learning a low-dimensional projection and the related reduced image from the noisy HSI. Then, the non-local low-rank denoising and iterative regularization are developed to refine the reduced image and projection, respectively. Finally, experiments on synthetic and both real datasets demonstrate the superiority against the other state-of-the-arts HSI denoising methods.READ FULL TEXT VIEW PDF
Matlab Code for: "Non-local Meets Global: An Integrated Paradigm for Hyperspectral Denoising. Arvix. Dec 2018"
Matlab Code for: "Non-local Meets Global: An Integrated Paradigm for Hyperspectral Denoising. Arvix. Dec 2018"
We apply tensor decoding scheme to decode SCI cameras.
Recent decades have witnessed the development of hyperspectral imaging techniques [4, 38, 18]. The hyperspectral imaging system is able to cover the wavelength region from 0.4 to 2.5 at a nominal spectral resolution of 10 nm. With the wealth of available spectral information, hyperspectral images (HSI) have the high spectral diagnosis ability to distinguish precise details even between the similar materials [2, 31], providing the potential advantages of application in remote sensing [32, 33], medical diagnosis 27, 33], quality control  and so on. Due to instrumental noise, HSI is often corrupted by Gaussian noise, which significantly influences the subsequent applications. As a preprocessing, HSI denoising is a fundamental step prior to HSI exploitation [6, 41, 43].
For HSI denoising, the spatial non-local similarity and global spectral correlation are the most two important properties. The spatial non-local similarity suggests that similar patches inside a HSI can be grouped and denoised together. The related methods [9, 12, 13, 28, 36, 45] denoise the HSIs via group matching of full band patches (FBPs, stacked by patches at the same location of HSI over all bands) and low-rank denoising of each non-local FBP group (NLFBPG). These methods have achieved state-of-the-art performance. However, they still face a crucial problem. For HSIs, the higher spectral dimension means the higher discriminant ability , thus more spectrums are desired. As the spectral number increases, the size of NLFBPG also becomes larger, leading to significantly more computations for the subsequent low-rank matrix/tensor approximations.
The HSIs have strong spectral correlation, which is modeled as low-rank property [1, 5, 41] and have also been widely adopted to the HSI denoising. However, only spectral low-rank regularization cannot remove the noise efficiently. One promising improvement is to project the original noisy HSI onto the low-dimensional spectral subspace, and denoise the projected HSI via spatial based methods [10, 29, 46, 47]. Unfortunately, these two-stage methods are significantly influenced by the quality of projection and the efficiency of spatial denoising. All of them fail to capture a clean projection matrix, which makes the restored HSI still be noisy.
To alleviate the aforementioned problems, this paper introduces a unified HSI denoising paradigm to integrate the spatial non-local similarity and global spectral low-rank property simultaneously. We start from the point that the HSI should underlies a low-dimensional spectral subspace, which has been widely accepted in hyperspectral imaging [16, 23], compressive sensing [3, 44], unmixing  and dimension reduction  tasks. Inspired by this fact, the whole NLFBPGs should also underlie a common low-dimensional spectral subspace. Thus, we first learn a global spectral low-rank projection, and subsequently exploit the spatial non-local similarity of projected HSI after the projection. The computational cost of non-local processing in our paradigm will almost keep the same with more spectral bands, and the global spectral low-rank property will also be enhanced. The contributions are summarized as follows:
We introduce a unified paradigm to exploit the spatial non-local and global spectral low-rank properties simultaneously. We transfer the non-local denoising to the reduced image and improve the computational efficiency against the increase of spectral band number;
The resulting new model for image denoising is hard to optimize, as it involves with both complex constraint and regularization. We further propose an efficient problem for optimization, which is inspired by alternating minimization;
Finally, the proposed method is not only the best compared with other state-of-the-art methods in simulated experiment, where Gaussian noise are added manually; but also achieves the most appealing recovered images for real datasets.
Notations We follow the tensor notation in , the tensor and matrix are represented as Euler script letters, i.e. and boldface capital letter, i.e. , respectively. For a -order tensor , the mode- unfolding operator is denoted as . We have , in which is the inverse operator of unfolding operator. The Frobenius norm of is defined by . The mode- product of a tensor and a matrix is defined as , where and .
Since denoising is an ill-posed problem, proper regulations based on the HSI prior knowledge is necessary [15, 35]. The mainstream of HSI denoising methods can be grouped into two categories: spatial non-local based methods and spectral low-rank based methods.
HSIs illustrate the strong spatial non-local similarity. After the non-local low-rank modeling was first introduced to HSI denoising in , the flowchart of the non-local based methods become fixed: FBPs grouping and low-rank tensor approximation. Almost all the researchers focused on the low-rank tensor modeling of NLFBPGs, such as tucker decomposition , sparsity regularized tucker decomposition , Laplacian scale mixture low-rank modeling , and weighted low-rank tensor recovery  to exploit the spatial non-local similarity and spectral low-rank property simultaneously. However, with the increase of spectral number, the computational burden also increases significantly, blocking the application of these methods to the real high-spectrum HSIs.
Chang et.al  claimed that the spectral low-rank property of NLFBPGs is weak and proposed a unidirectional low-rank tensor recovery to explore the non-local similarity. It saved much computational burden and achieved the state-of-the-art performance in the HSI denoising. This reflects the fact that previous non-local low-rank based methods have not yet efficiently utilized the spectral low-rank property. How to balance the importance between spectral low-rank and spatial non-local similarity still remains a problem.
, the intrinsic dimension of the spectral subspace is far less than the spectral dimension of the original image. By vectorizing each band of the HSI and reshaping the original 3-D HSI into a 2-D matrix, various low-rank approximation methods such as principal components analysis (PCA), robust PCA [11, 37, 41], low-rank matrix factorization [3, 39] have been directly adopted to denoise the HSI. However, these methods only explore the spectral prior of the HSI, ignoring the spatial prior information. Instantly, many conventional spatial regularizers such as total variation , low-rank tensor regularization [25, 30] are adopted to explore the spatial prior of HSI combined with spectral low-rank property.
A remedy is a two-stage method combining the spatial regularizer and spectral low-rank property together. This is done by firstly mapping the original HSI into the low-dimensional spectral subspace, and then denoise the mapped image via existing spatial denoising methods, e.g., wavelets [10, 29], BM3D  and HOSVD . These two-stage methods provide a new sight to denoise the HSI in the transferred spectral space, which is very fast. However, these methods fail to combine the best of both worlds, and the extracted subspace is still corrupted by the noise.
In this section, we propose a unified HSI denoising paradigm to integrate spatial non-local similarity and global spectral low-rank property. We first learn a low-dimensional projection and the related reduced image from the noisy HSI. Then the reduced image and the projection are updated by spatial non-local denoising and iteration regularization, respectively. The overview of the proposed paradigm is in Figure 1.
Assuming that the clean HSI is corrupted by the additive Gaussian noise
(with zero mean and variance), then the noisy HSI is generated by
First, to capture the spectral low-rank property in Section 2.2, we are motivated to use a low-rank representation of the clean HSI , i.e. , where , is a projection matrix capturing the common subspace of different spectrum, and is the reduced image. Second, to utilize the spatial low-rank property, we add a non-local low-rank regularizer on the reduced image . As a result, the proposed non-local meets global (NGmeet) denoising paradigm is presented as
where controls the contribution of spatial non-local regularization, the projection matrix is required to be orthogonal, and the clean HSI is recovered by .
The objective (2) is very hard to optimize, due to both the orthogonal constraint on and complex regularization on . An algorithm based on alternating minimization to approximately solve the objective function is proposed in Section 3.2.
The orthogonal constraint is very important here.
it encourages the representation held
in to be more distinguish with each other.
This helps to keep noise out of and
further allows a closed-form solution for computing (Section 3.2.1 ).
it preserves the distribution of noise,
which allows us to estimate the remained noise-level in reduced image and
reuse the-state-of-art Gaussian based non-local method for spatial denoising (Section
). Besides, it preserves the distribution of noise, which allows us to estimate the remained noise-level in reduced image and reuse the-state-of-art Gaussian based non-local method for spatial denoising (Section3.2.2).
Recall that, in (2), the first item tries to exploit the spectral low-rank property and decompose the noisy into the coarse spectral low-rank projection and reduced image . Here, both and have physical meaning in the field of remote sensing . Specifically, -th column of , denoted as , is regarded as the -th signature (known as endmember) of HSI, and the corresponding coefficient image is regarded as the abundance map.
Previous methods are mostly two-stage ones, they do not iterative refine the projection matrix they found, e.g. FastHyDe . However, we model the spatial and spectral low-rank properties simultaneously, which enables iterative refinement of the projection matrix . To demonstrate the necessity of iterative refinement, we calculated the projection and reduced image from noisy WDC with noise variance 50. The reference and are from the original clean WDC. Figure 2 presents the comparison on signatures and the corresponding coefficient image before and after our refinement. From the figure, it can be observed that the projection atom and reduced image obtained by the spectral denoising method are still suffering from the noise, while the proposed method can produce much cleaner signatures and coefficient images.
As discussed in Section 3.1, the objective (2) is very hard to optimize. In this section, we are motivated to use alternating minimization for optimization (Algorithm 1). , stand for the input noisy image and output denoised image of the -th iteration, respectively. As will be shown in the sequel, Algorithm 1 tries to find a closed-form solution for (step 3) and reuses the-state-of-art spatial denosing method for computing (steps 4-6), which together make the algorithm very efficient. Besides, as will be refined during the iteration, iterative regularization  is adopt to boost the denosing performance (step 7).
In this stage, we identify the projection matrix with the given and from (2), which leads to
However, this problem is hard without simple closed-form solution. Instead, since is obtained from iterative regularization, of which the noisy-level is decreased. Thus, we proposed to relax (3) as
which has the closed-form solution (Proposition 3.1). Thus, only a SVD on the folding matrix of is required, which can be efficiently computed.
Let be the rank- SVD of . The solution to (4) is given by the close-form as and .
where is a non-local denoising regularizer. Formulation (5) appears in many denoising models, e.g. WNNM . Specifically, to solve this regularizer, we need to first group similar patches, then denoise each patch group tensors and finally assemble the final estimated .
However, all these model assume the noise on
follow univariate Gaussian distribution. If such assumption fails, the resulting performance can deteriorate significantly. Here, we have the following Proposition3.2. Therefore, the noise distribution is preserved from to , which enables us to use the existing spatial denoising methods. In this paper, we use WNNM  to denoise each patch group tensor, as it is widely used and gives the-state-of-art denoising performance.
Assume the noisy HSI is from (1), then the noise on the reduced image , where , still follows Gaussian distribution with zero mean and variance .
Finally, to use WNNM, we need to estimate the noise level in , whose noise level is changed during the iteration. From Proposition 3.2, we know the noisy level of is the same as , thus we propose to estimate it via
where is the a scaling factor controlling the re-estimation of noise variance, and stands for the averaging process of the tensor elements. The denoised group tensors are denoted as , which can be directly used to reconstruct the denoised reduced image . The output denoised image of -th iteration is .
Iteration regularization has been widely used to boost the denoising performance [9, 14, 19, 36]. Here we also introduce the it into our model (Algorithm 1) to refine the noisy projection . As shown in (4), the projection is significantly influenced by the noise intensity of input noisy image . Hence we update the next input noisy image as
where is to trade-off the denoised image and original noisy image . The estimation of can benefit from the lower noise variance of the input .
Besides, is also updated with the iteration. We initialize by HySime . When the noisy image is corrupted by heavy noise, the estimated
will be small. Fortunately, the larger singular values obtained from the noisy image are less contaminated by the noise, and help to keep noise out of the reduced image. With the iteration, We increaseby
where is a constant value. Therefore, has the ability to capture more useful information with more iterations.
|stage A||stage B|
Following the procedure of Algorithm 1, the main time complexity of each iteration includes stage A-SVD , stage B.non-local low-rank denoising of each . Table 1 presents the time complexity comparison between NGmeet and other non-local HSI denoising method. LLRT and KBR only need stage B to complete the denoising. As can be seen, the proposed NGmeet costs additional complexity in stage A, however, will be at least times faster in stage B.
|spectral low-rank methods||spatial low-rank methods|
In this section, we present the simulated and real data experimental results of different methods, companied with the computational efficiency and parameter analysis of the proposed NGmeet. The experiments are programmed in Matlab with CPU Core i7-7820HK 64G memory.
Setup. One multi-spectral image (MSI) CAVE 111http://www1.cs.columbia.edu/CAVE/databases/, and two HSI images, i.e. PaC 222http://www.ehu.eus/ccwintco/index.php/ and WDC 333https://engineering.purdue.edu/~biehl/MultiSpec/hyperspectral datasets are used (Table 3). These images have been widely used for a simulated study [9, 21, 28, 36, 47]. Following the settings in [9, 28], additive Gaussion noise with noise variance are added to the MSIs/HSIs with varies from to . Before denoising, the whole HSIs are normalized to [0, 255].
|number of bands||31||89||192|
The following methods are used for the comparison: spectral low-rank methods, i.e. LRTA  444https://www.sandia.gov/tgkolda/TensorToolbox/, LRTV  555https://sites.google.com/site/rshewei/home, MTSNMF  666http://www.cs.zju.edu.cn/people/qianyt/, NAILRMA  PARAFAC  and FastHyDe  777http://www.lx.it.pt/~bioucas/; spatial low-rank methods, i.e. TDL  KBR  888http://gr.xjtu.edu.cn/web/dymeng/, LLRT  999http://www.escience.cn/people/changyi/; and finally NGmeet (Algorithm 1), which combines the best of above two fields. Hyper-parameters of all compared methods are set based on authors’ codes or suggestions in the paper. The value of spectral dimension is the most import parameter, which is initialized by HySime  and updated via (7). Parameter is used to control the contribution of non-local regularization, and is a scaling factor controlling the re-estimation of noise variance . We empirically set , and as introduced in , and in the whole experiments.
To thoroughly evaluate the performance of different methods, the peak signal-to-noise ratio (PSNR) index, the structural similarity (SSIM)  index and the spectral angle mean (SAM) [9, 22] index were adopted to give a quantitative assessment. The SAM index is to measure the mean spectrum degree between the original HSI and the restored HSI. The lower value of SAM means the higher similarity between original image and the denoised image.
Quantitative comparison. For each noise level setting, we calculate evaluation values of all the images from each dataset, as presented in Table 2. It can be easily observed that the proposed NGmeet method achieved the best results almost in all cases. Another interesting observation is that the non-local based method LLRT can achieve better results than FastHyDe, the best result of spectral low-rank methods, but it dose the opposite in the hyperspectral image cases. This phenomenon conforms the advantage of NL low-rank property in the MSI processing and the spectral low-rank property in the HSI processing.
Visual comparison. To further demonstrate the efficiency of the proposed method, Figure 3 shows the color images of CAVE-toy (composed of bands 31, 11 and 6 ) before and after denoising. The results of PaC and WDC can be found in the supplementary material. The PSNR values and the computational time of each methods are marked under the denoised images. It can be observed that FastHyDe, LLRT and NGmeet have huge advantage over the rest comparison methods. From the enlarged area, the results of FastHyDe LLRT produced some artifacts. Thus, our method NGmeet can produce the best visual quality.
|(seconds)||stage B||stage B||stage A||stage B||total|
Computational efficiency. In this section, we will illustrate that in our denoising paradigm, the computational efficiency of the non-local denoising procedure will get rid of the huge spectral dimension. Compared to the previous non-local denoising methods, i.e. KBR  and LLRT , the proposed NGmeet includes additional stage A. Table 4 presents the computational time of different stages of the three methods. From Table 1 and 4, we can conclude that NGmeet spends little time to project the original HSI onto a reduced image (stage A), however, earning huge advantage in stage B including group matching step and non-local denoising.
Figure 4 displays the computational time and SSIM values of the proposed NGmeet, KBR  and LLRT , with the increase of spectral number. As illustrated, even though the performances of KBR and LLRT increase with the increase of spectral number, the computational time also increases linearly. Our method can achieve the best performance, meanwhile, the computational time is nearly unchanged with the increase of spectral number.
Convergence. To show the convergence of proposed NGmeet, Figure 5 presents the PSNR values with the increase of iteration, on the WDC dataset. From the figure, it can be observed that our method can converge to stable PSNR values very fast at different noise level.
Setup. Here, AVIRIS Indian Pines HSI 101010https://engineering.purdue.edu/~biehl/MultiSpec/ and HYDICE Urban image 111111http://www.tec.army.mil/hypercube are adopted in the real experiments (Table 5). The noisy HSIs are also scaled to the range [0 255], and the parameters involved in the proposed methods are set as the same in the simulated experiments. In addition, multiple regression theory-based approach  is adopted to estimate the initial noise variance of each HSI bands.
|number of bands||210||220|
Visual comparison. Since reference clean images are missing, we just present the real Indian Pines and Urban images before and after denoising in figures 6 and 7. It can be obviously observed that the results produced by the proposed NGmeet can remove the noise and keep the spectral details simultaneously. LRTV can produce the most smooth results. However, the color of the denoised result changes a lot, indicating the loss of spectral information. The denoised results of FastHyDe and LLRT still contain stripes as presented in Figure 6. To sum up, although the proposed NGmeet is designed in the Gaussion noise assumption, it can also achieves the best results for real datasets.
is the key parameter to integrate the spatial and spectral information. Figure 8 presents the PSNR values achieved by NGmeet with different initialization of with being . PaC images was chosen as the test image, and the noise variance changes from to . is initialized by HySime  as for different noise variance cases, respectively. It confirms that the initialization of is reliable.
In this paper, we provide a new perspective to integrate the spatial non-local similarity and global spectral low-rank property, which are explored by low-dimensional projection and reduced image denoising, respectively. We have also proposed an alternating minimization method with iteration strategy to solve the optimization of the proposed GNmeet method. The superiority of our method are confirmed by the simulated and real dataset experiments. In our unified spatial-spectral paradigm, the usage of WNNM 
is not a must. In future, we plan to adopt Convolutional Neural Network[7, 42] to explore non-local similarity.
Hyperspectral image super-resolution via non-local sparse tensor factorization.In , pages 3862–3871, July 2017.
Note that the objective can be expressed as:
which is equal to find the best -rand approximation of . Thus, let rank- SVD of be , and , the closed-form solution of (4) is given by and . ∎
Since , then
where the noise is given by . Note that
Thus, the mean of the noise is zero. Let be a column in , then one column in can be expressed as
Follow the definition of variance, we have
Thus, we obtain the proposition. ∎
Very recently deep learning based HSI denoising mehtods,i.e. HSID-CNN  and HSI-DeNet  have been proposed. We didn’t compare the proposed NGmeet with these two methods, since the related codes are no public available. However, from the cross-over experimental report between our method and the two deep learning based methods, we found that our method performs much better than HSID-CNN  in the WDC data case with Gaussion noise. HSI-DeNet  didn’t report the Gaussion noise removal of CAVE dataset. However, if we regard LLRT as a baseline, the inprovement of our NGmeet compared to LLRT is a bit higher than that of HSI-DeNet . These evidences prove the validity of our method compared to the state-of-the-art deep learning based methods.