Fusion of Hyperspectral and Multispectral Images Based on a Centralized Non-local Sparsity Model of Abundance Maps

INTRODUCTION

Hyperspectral (HS) imaging acquires a real scene in hundreds of spectral bands, each associated with a specific wavelength (Medina et al., 2017). The analysis of spectral signatures of these images has allowed the advancement of many applications in the fields of medical imaging (Panasyuk et al., 2007; Lu and Fei, 2014), remote sensing (Schaepman et al., 2009; Bioucas-Dias et al., 2013; Velasco et al., 2016), and astronomy (Schreiber et al., 2009; Hege et al., 2004). However, due to hardware restrictions, HS images are limited by a low spatial resolution (Shaw and Burke, 2003). For example, the Hyperion imaging spectrometer has a spatial resolution of 30 meters-per-pixel (Middleton et al., 2013) that can degrade the performance for practical applications.

In order to increase the spatial resolution of HS images, a usual approach is to merge the HS image with high spatial resolution images. An example appropriate for this paper consists of merging the HS image (that has high spectral resolution) with a multispectral (MS) image (that has high spatial resolution) (Wei et al., 2015; Simoes et al., 2015). Another well-known example is HS pan-sharpening, which addresses the fusion of panchromatic and HS images (Loncan et al., 2015).

Image fusion is essential to find a suitable model considering the prior knowledge of natural images due to the ill-posed nature of the inverse problem. Thus, some effective regularizers that restrict the solution space have been employed in the state-of-the-art for image restoration (IR) problems with encouraging results (Dong, Zhang and Shi, 2011; Chan et al., 2015; Oliveira, Bioucas-Dias and Figueiredo, 2009; Vargas et al., 2018). In detail, under a super-resolution approach, employing previous total variations have achieved state-of-the-art performance for fusing HS and MS images (Simoes et al., 2015). Additionally, it has been demonstrated that assuming the sparcity of the signal of interest is a suitable regularization method in image restoration problems such as image de-blurring, image denoising, and image super-resolution (Dong et al. 2011, Wei et al., 2015; Yang et al., 2010). Formally, a sparce representation states that a signal 𝐱 ∈ ℝ^Np can be represented as a linear combination of a few atoms via 𝐱 = 𝚽𝜶_x on the dictionary Φ (Cardon et al., 2015). If the target signals reconstructed are multi-dimensional images, the representation of smaller patches on an over-complete dictionary learned from the data successfully determines the underlying structure of the images (Mallat, 1999).

On the other hand, restoring the image x from a degraded image y = Hx + n under sparcity is challenging due to degradation during the acquisition process. The estimation of the sparce coefficient cl _x employing common optimization algorithms as Ba-yesian frameworks, iterative shrinkage/thresholding, Lasso, and basis pursuit (Tibshirani, 1996; Tropp and Wright, 2010; Beck and Teboulle, 2009; Bruckstein, Donoho and Elad, 2009; Chen, Donoho and Saun-ders, 2001) may not be an optimal estimation of the original sparce coefficient a _x . Thus, the reconstruction may lead to an inaccurate estimation (Mairal et al., 2009; Dong et al., 2012).

According to this limitation, we are interested on improving the estimated sparce representation, so non-local redundancy properties in natural scenes are studied to enhance the representation model. For instance, non-local self-similarities of natural scenes have been included in approaches for super-resolution or denoising leading to state-of-the-art performance (Elad and Aharon; 2006; Dian, Fang and Li, 2017; Mairal et al., 2009; Fu et al., 2017). Similarly, to enhance the accuracy of IR methods based on sparsity, a centralized non-local sparce representation (CNSR) model has been proposed in (Dong et al., 2012). Specifically, this representation intends to increase the quality of the reconstructed image by reducing the error v _α = α _x - α _y known as sparce coding noise (SCN). This goal is achieved by centralizing the sparce coefficient α _y to α _x that is estimated from the non-local redundancies in x.

On the other hand, research on IR literature has included the unmixing spectral model on several IR problems to find a treatable connection between the distorted measurements and the original image, and to take advantage of intrinsic properties such as spectral non-local similarity (Wei et al., 2016; Vargas, Arguello and Tourneret, 2019; Zhao, Yang and Chan, 2013). Even sparcity assumptions in abundance (Vargas et al., 2018, Ghasrodashti et al., 2017; Zhao et al., 2014) and spatial non-local similarity (Zhao, Yang and Chan, 2013) have been taken into account on the hyperspectral image super-resolution and fusion problems showing promising results. Nonetheless, prior unified knowledge of spatial-spectral non-local similarity and sparcity has not been considered on the abundances.

This work introduces a multispectral and hyperspectral image fusion model that assumes a CNSR of abundance maps in a learned dictionary. More precisely, local and spectral non-local similarity assumptions on the abundances are included in the model by the spectral unmixing decomposition and the spatial non-local estimation from α _y now associated with the sparse representation of the abundances α _α , respectively. Additionally, the abundances are supposed to be sparce on an over-complete dictionary, which is constructed from available abundance and HS-MS data. The dictionary consists of several sub-dictionaries learned from groups of similar patches. Then, all patches of a cluster are sparsely represented and centralized to a good estimation of α _α under its corresponding sub-dictionary (each cluster associates with a single sub-dictionary). Therefore, both the construction of a composite over-complete dictionary and the centralization are performed to exploit non-local redundancies in the abundances.

In the global optimization problem, the sparse codes and spatial-spectral dictionary are estimated jointly under an alternating approach. Once the dictionary has been updated, an iterative nonlinear shrinkage algorithm solves the resulting multi-band image fusion problem; this algorithm consists of projecting the solution of a quadratic problem on the thresholding operator in each iteration (Daubechies, Defrise and De Mol, 2004). Moreover, the regularization parameters underlying the optimization problem are adaptively adjusted from a Bayesian formulation (Dong et al. 2012). Non-local similarity and sparcity assumptions are favorable on abundance maps due to the large presence of soft regions. Additionally, the burden load decrease by data reduction compared to the original image. Considering experimental scenarios with real HS images, we show that the proposed fusion approach outperforms and achieves the performance of the competitive approaches gaining up to 1[dB] under the peak signal-to-noise ratio (PSNR) metric.

PROBLEM STATEMENT

For modeling the MS and HS images, it is commonly assumed that there is a linear degradation in the spectral and spatial domain of a HS image with higher resolution denoted by the vector with N _p =NML (Wei et al., 2016; Wei, Dobigeon and Tourneret, 2015; Simoes et al., 2015). The acquired HS image denoted by the vector is assumed to be a blurred and down-sampled version of the high full-resolution image x. In contrast, the MS image is obtained from a spectral degradation of the target image (Wei et al., 2015). Thus, the sensing models for the MS and HS images are given as

Where represents the down-sampling matrix, is a circular convolution operator described as blurring matrix, denotes the transfer function of the MS sensor, and the vectors and represent additive Gaussian noise for the MS and HS images.

For notation convenience, the spatial-spectral tensor (L is the number of spectral bands and N x M represents the spatial dimension) is represented by the vector , where contains all the image intensities associated with the z'-th spectral band on X. Likewise, we can exploit intrinsic properties of the spectral unmixing model for the fusion problem in which each spectral signature of a multi-band image can be associated to a linear mixture of several pure spectral signatures known as endmembers. Under this model, the image x is decomposed as ^x =Ma, where with N _L =NML are the P abundane maps in a vector form denoting as the z'-th abundance map. The matrix is a block-diagonal matrix on the endmembers' matrix with P the number of materials. Thus, the acquisition models in equation (1) can be rewritten as

Considering the observation models in equation (2), the optimization problem formulated in this paper consists of estimating the vector of abundances from the observed measurements y _M and y _H . Since the resultant problem is ill-posed, this work proposes that the regularization strategy employed for the fused solution be a centralized non-local sparce model on the abundance maps.

METHODOLOGY

In this section, we describe the methodology for the fusion of hyperspectral and multispectral images based on the non-local sparse representation model. First of all, consider that the abundance data a is decomposed into overlapping image cubic patches. Dividing images in patches has obtained effective results in the image domain in many image processing applications (Tosic and Frossard, 2011). We start by defining as the matrix extracting cubic patch at location as the vectorized abundance cubic patch with dimensions S x S x P, and i indicating the position of the central pixel. Then, considering a sparce representation of the cubic patch in a fixed dictionary, i.e., = Φ α _α,i , where α _α,i is a sparce vector. Finally, considering the collection of sparce coefficients {α _α,i }, the recovery of the abundances a can be computed by averaging all the patches

For notation simplicity, the estimated abundance ậ aforementioned is denoted as ậ =Φ o α _α , where α _α represents the collection of the set of sparce codes α _α,i .

Dictionary learning

Basis functions extracted from learned or analytical dictionaries are able to represent basic structures of natural image and define sparcity domains. Let a given cubic patch a _i be represented as a matrix where the rows are the abundance coefficients. Under a sparce model of the abundances, A _i , can be described by few atoms (Fu et al., 2017)

where {Φ_m} denotes a collection of atoms, α _m are the corresponding sparce codes, and m « S ² P. The selection of the dictionary elements is a relevant issue of sparcity-based methods. In this work, we propose a dictionary with tridimensional atoms that are constructed using separable abundance and spatial components as follows

where are orthonormal basis spanning the abundance maps and spatial subs-pace, respectively. Under this factorization of the dictionary, we are able to exploit the correlation in the abundances as well as in the spatial domain. Specifically, the basis for the abundance subspace domain U and the spatial domain V are extracted from the abundances data allowing to construct the atoms in equation (12) with correlation information between abundance maps and spatial information simultaneously. In the abundance subspace domain, the estimated abundance vector â isreorganized as a P x NM matrix via PCA the abundance basis U = [u ₁ , -, Up]. Similarly, the spatial basis V = [v ₁ , …, v_S2] is calculated using PCA from spatial patches extracted along all the abundance maps.

In the spatial domain, we learn different dictionaries for different clusters of the abundance training patches (Dong et al., 2011 ; Dong et al., 2012). Thus, the patches are clustered by spatial similarity of â pooled along with all the abundance maps employing the K-means method. This methodology lies on grouping patches spatially in different clusters C _K , with k = 1, … ,K, by using the high frequency components of the abundance maps. This strategy is more suitable since clustering by intensity presents low performance discriminating images (Fu et al., 2017). The high frequency patterns ã of â can be calculated as

where G is a low pass filtering operator. It is worth noting that the high frequency of the abundance maps ã is employed only for clustering and the resultant clusters C _K associates with a spatial basis extracted from the estimated abundance maps â. Composing the PCA abundance and spatial basis following equation (12), the proposed structured spatial-spectral dictionary for a given abundance cubic patch a _i in C _K can be written as follow

where cienotes the Kronecker operator. It is important to note that a _i only can be associated with a single spatial-spectral dictionary . A schematic representation of the spatial-abundance structure of the proposed adaptive dictionary is shown in Figure 1.

Figure 1: Schematic representation of the proposed adaptive spatial-spectral dictionary learned from the abundance maps.

Source: Authors.

Fusion algorithm based on a CNSR model

The strategy to solve the MS and HS image fusion problem based on a non-local sparce representation of the abundance maps in equation (4) is an iterative scheme that employs a shrinkage step in every iteration (Dong et al., 2012). The proposed algorithm is summarized in Algorithm 1.

Algorithm 1:: MS and HS Image Fusion based on CNSR

The algorithm consists of two main iterative loops. The external loop (whose iterations are indexed by t) is implemented to update the spatial-abundance dictionary . Once the dictionary is computed, the internal loop (whose iterations are indexed by q) estimates an approximation of the abundances â by employing an alternating strategy of three steps. First, the solution of the problem associated with the forward model or the image fidelity term is found by a gradient descent algorithm. Second, a rough estimation of the abundance maps is obtained by the sparce unmixing algorithm via variable splitting and augmented Lagrangian (SUnSAL) introduced in (Bioucas et al., 2010). Third, a shrinkage operator (Dong et al., 2012; Daubechies, Defrise and De Mol, 2004) is employed to compute the sparce coefficients â_y . Thus, the solution â_y in the q + 1-th iteration is given by

where is the soft-thresholding proximal solution of the ℓ1 norm with a threshold parameter t.

EXPERIMENTAL RESULTS

This section presents the performance of the proposed method on two target hyperspectral images fixed as ground truth for the simulations. The experimental scenarios are conducted on both data sets: the Pavia image with spatial resolution 128 x 128 and 93 spectral bands, and the Moffett image with spatial resolution 128 x 128 and 93 spectral bands. The high-resolution ground-truth images are degraded by spectral and spatial distortion operators to obtain the MS and HS images, respectively. Particularly, the HS image is obtained from the target image by filtering every band with a 5 x 5 filter sampled from a Gaussian distribution followed by a down-sampling operation on the vertical and horizontal directions with a scale factor d = 4. The MS image was calculated by applying a spectral filtering of the high-resolution ground-truth images with the LANDSAT g spectral response leading to a MS image with L _M = 4 bands. Then, both degraded MS and HS images are perturbed band-by-band with additive Gaussian noise. The HS image is contaminated with SNR 30 = [db] in the last 50 bands and SNR 35 = [db] for the bands remaining while for the MS image is contaminated with a SNR 30 = [db]. The setting parameters for the fusion algorithm were selected as follows: K = 60 clusters; spatial patch size 6 x 6, 12 similar neighbors cubic patches are extracted to calculate β _i (i.e., |Ωi|) in a neighborhood size of 50 x 50; the parameter δ is set on 1.25 while the thresholding for sparce representation is set on t = 1.841. In order to compare the proposed method, three different approaches known as (HMIF-SR), (HySure), and (HyMuXNol) (Wei et al., 2015; Simoes et al., 2015; Arias, Vargas and Arguello, 2019) are simulated on the same fusion scenario.

The MS and HS images and the full-resolution reconstructed images for the Pavia and Moffett images are shown in Figures 2 and 3. The results from the proposed fusion strategy show a better-quality estimation of the target image. Furthermore, the fused resultant images of adversary methods depict greater spectral and spatial degradation compared to the Pavia and Moffett ground-truth images.

Figure 2: Qualitative fused images (Pavia image). (Row 1, left) MS image. (Row 1, middle) MS image. (Row 1, right) Ground truth. (Row 2, left) HySure. (Row 2, middle) HMIF-SR. (Row 2, right) HyMuXNol. (Row 3) Proposed.

Source: Authors.

Figure 3: Qualitative fused images (Moffet image). (Row 1, left) MS image. (Row 1, middle) MS image. (Row 1, right) Ground truth. (Row 2, left) HySure. (Row 2, middle) HMIF-SR. (Row 2, right) HyMuXNol. (Row 3) Proposed.

Source: Authors.

The evaluation of the proposed fusion method is performed employing the PSNR (Peak Signal-to-Noise Ratio), UIQI (Universal Image Quality Index), RMSE (Root-Mean-Square error), ERGAS (Relative Dimensionless Global Error in Synthesis), SAM (Spectral Angle Mapper), and DD (Degree of Distortion) image fusion metrics (see (Wei et al., 2015) for more details). Tables 1 and 2 summarize the computed quantitative metrics of the proposed strategy compared to the competitive fusion strategies whose numerical values clearly favor the restored images from the suggested fusion method.

Table 1: Quantitative fusion results of the proposed method and the competitive algorithms on the Pavia image: PSNR [dB], UIQI, RMSE (102), ERGAS, SAM [degrees], DD (103)

Source: Authors.

Table 2: Quantitative fusion results of the proposed method and the competitive algorithms on the Moffet image: PSNR [dB], UIQI, RMSE (102), ERGAS, SAM [degrees], DD (103)

Source: Authors

To further show the advantage of the proposed solution compared to the best method in the state-of-the-art HyMuXNol, we randomly extracted four spectral signatures on the Pavia and Moffett ground-truth images and evaluated them qualitatively (Figure 4). It is possible to see that the recovered spectral signatures using the prior CNSR on the abundances are closer to the ground-truth spectral signatures, indicating that the spectral information of the scenes is better fitted with the proposed model than the underlying model in the HyMuXNol method, and this happens for all the selected images.

Figure 4: Evaluation of the spectral performance: (Row 1) Ground-truth and recovered spectral signatures on the Pavia data set. (Row 2) Ground-truth and recovered spectral signatures on the Moffett data set.

Source: Authors.

CONCLUSIONS

This work proposed a MS and HS image fusion model solved under a sparce model on a set of learned dictionaries. By using the spectral unmixing decomposition, this model includes a non-local sparce representation of abundance maps on the fusion problem while taking advantage at the same time of the low dimensionality of the HS data. The estimated sparce representation is centralized for a better estimation obtained from the non-local redundancy of the abundances. These considerations lead to enhance the quality of fused images of state-of-the-art fusion methods based on sparce regularization. Additionally, a spatial-spectral dictionary is adaptively constructed, exploiting the low dimensionality of the abundance maps. The dictionary structure consists of K sub-dictionaries, each one estimated from a cluster of similar spatial features of abundance maps. Thus, an abundance patch is then sparcely represented in the appropriate dictionary. In order to obtain the full resolution image, the employed numerical strategy includes two main steps alternated iteratively: shrinkage thresholding operator for sparse regularization is used to solve the ℓ1 norm, and the execution of a gradient descent method to solve the quadratic fidelity problem. Restoration results show that the proposed image fusion model based on the abundance map analysis outperforms the competitive fusion methods based on sparse regularization in terms of the quantitative metrics PSNR, UIQI, RMSE, ERGAS, DD, and SAM.