# Hyperspectral Pansharpening Based on Homomorphic Filtering and Weighted Tensor Matrix

^{1}

^{2}

^{3}

^{*}

Joint Laboratory of High Speed Multi-Source Image Coding and Processing, School of Telecommunications Engineering, Xidian University, Xi’an 710071, China

State Key Lab. of Integrated Service Networks, School of Telecommunications Engineering, Xidian University, Xi’an 710071, China

Department of Electrical and Computer Engineering, Mississippi State University, Mississippi State, MS 39762, USA

Author to whom correspondence should be addressed.

Received: 5 February 2019 / Revised: 22 April 2019 / Accepted: 23 April 2019 / Published: 27 April 2019

Hyperspectral pansharpening is an effective technique to obtain a high spatial resolution hyperspectral (HS) image. In this paper, a new hyperspectral pansharpening algorithm based on homomorphic filtering and weighted tensor matrix (HFWT) is proposed. In the proposed HFWT method, open-closing morphological operation is utilized to remove the noise of the HS image, and homomorphic filtering is introduced to extract the spatial details of each band in the denoised HS image. More importantly, a weighted root mean squared error-based method is proposed to obtain the total spatial information of the HS image, and an optimized weighted tensor matrix based strategy is presented to integrate spatial information of the HS image with spatial information of the panchromatic (PAN) image. With the appropriate integrated spatial details injection, the fused HS image is generated by constructing the suitable gain matrix. Experimental results over both simulated and real datasets demonstrate that the proposed HFWT method effectively generates the fused HS image with high spatial resolution while maintaining the spectral information of the original low spatial resolution HS image.

Depending on the number of acquired bands, remote sensing imaging technology has developed from collecting panchromatic (PAN) and color images to multispectral (MS) images, and it can now capture hyperspectral (HS) images with dozens of hundreds of bands. A PAN image with very high spatial resolution is a single-band grayscale image acquired in the visible range. It is able to obtain the shape feature of objects, but cannot distinguish colors. A color image consists of three bands which are red, green and blue, and displays the colors of objects. However, it is difficult to distinguish the features in similar colors. An MS image not only obtains spatial features, but also obtains spectral information in several bands, which is more capable of distinguishing categories of different features. However, the rough spectral resolution of MS images may not meet the requirements in some applications, and it is hard to realize fine feature detection [1]. An HS image with a higher spectral resolution on the order of nanometers can provide finer classification [2], which has been applied to many fields [3,4,5,6,7] and some practical applications, such as vegetation study [8], precision agriculture [8], regional geological mapping [9], mineral exploration [10], and environment monitoring [11]. Due to technical limitations, the spatial resolution of an HS image is low.

As both high spatial and spectral resolutions are important in practical applications, obtaining a high spatial resolution HS (HRHS) image is crucial. One effective way is to perform hyperspectral pansharpening, which fuses a high spatial resolution PAN (HRPAN) image with a low spatial resolution HS (LRHS) image. Figure 1 shows the concept of hyperspectral pansharpening.

Many hyperspectral pansharpening algorithms were developed, among which hyperspectral pansharpening methods using Bayesian and matrix factorization have been proposed in recent years. The Bayesian-based approaches include Bayesian naive Gaussian prior [12], Bayesian sparsity promoted Gaussian prior [13], and HySure [14]. These algorithms utilize the posterior distribution, and are based on maximum a posteriori estimation to fuse LRHS and HRPAN images [15]. The matrix factorization approach generates a fused HRHS image by using the nonnegative matrix factorization (NMF) under some constraints to estimate endmember and abundance matrices [16]. The matrix factorization approach is well represented by the nonnegative sparse coding (NNSC) [17] and constrained nonnegative matrix factorization (CNMF) [18] methods. The main challenge in hyperspectral pansharpening is to effectively improve the spatial resolution while preserving the original spectral information. The Bayesian and matrix factorization approaches are able to achieve good results on this challenge, but have a high computational cost.

Component substitution (CS) and multi-resolution analysis (MRA) approaches are two classical hyperspectral pansharpening approaches which have simple and fast implementation. For the CS class, intensity-hue-saturation (IHS) transform [19,20], principal component analysis (PCA) transform [21,22], Gram–Schmidt (GS) [23], and adaptive GS (GSA) [24] are the most representative methods. The CS class extracts spatial details of the HS image, and replaces the extracted spatial details with the HRPAN image. Regardless of superior spatial performance, the CS class suffers from serious spectral distortion [25]. The typical algorithms of the MRA technique are smoothing filter based intensity modulation (SFIM) [26], Laplacian pyramid [27], modulation transfer function generalized Laplacian pyramid (MTF-GLP) [28], and MTF-GLP with high pass modulation (MTF-GLP-HPM) [29]. The MRA methods generally utilize a multi-resolution decomposition to extract spatial details which are imported into the HS image. Compared with the CS methods, the MRA methods generate less spectral distortion, but usually have a larger computational burden [30]. Recently, several algorithms based on the CS and MRA approaches have been proposed, such as the Sentinel-2A CS and MRA based sharpening algorithm [31], the multiband Filter estimation (MBFE) algorithm [32], and the guided filter PCA (GFPCA) algorithm [33]. Moreover, several intelligent processing-based methods have also been proposed, and examples include deep two-branches convolutional neural network (Two-CNN-Fu) [34], Bidirectional Pyramid Network [35], and 3D-convolutional neural network (3D-CNN) [36].

The CS and MRA approaches mostly extract the spatial information of the HRPAN image and inject it into the LRHS image, but without considering the spatial information of the LRHS image. Due to the incomplete spatial information injection, the CS and MRA approaches may result in distortion. To address this problem, we propose a novel hyperspectral pansharpening method by combining homomorphic filtering with a weighted tensor matrix. An optimized weighted tensor matrix-based method which considers the structure information of the LRHS and HRPAN images is proposed to generate more comprehensive spatial information. In addition, to extract the spatial structure information of the LRHS images, open-closing morphological operation is first used for noise removal, and homomorphic filtering is then introduced to extract the spatial details of each band. Finally, a weighted root mean squared error based method is proposed to obtain the total spatial component of the LRHS image from extracted spatial details of each band, and the Laplacian pyramid networks super-resolution algorithm is adopted to enhance the spatial resolution of the obtained spatial component. Comparative analysis was used to demonstrate the applicability and superiority of the proposed method in both spectral and spatial qualities.

As stated above, a new hyperspectral pansharpening method based on homomorphic filtering and weighted tensor matrix is proposed in this paper. The main novelties of the proposed hyperspectral pansharpening method are concluded in the following aspects.

- A novel HS image spatial component extraction strategy is proposed. Open-closing morphological operation and homomorphic filtering are first introduced to remove the noise and extract the spatial details of each band of the HS image, respectively. Then, a weighted root mean squared error-based method is proposed to obtain the total spatial component of the HS image.
- An optimized weighted tensor matrix-based method is proposed to integrate the spatial component of the HS image with the spatial component of the PAN image. The weighted structure tensor matrix that represents the structural information of multiple images is applied to hyperspectral pansharpening for the first time. The classical methods which mostly extract the spatial information of the PAN image inject the incomplete spatial information, and may lead to distortion. Unlike there classical methods, the proposed optimized weighted tensor matrix-based method generates the spatial information not only from the PAN image but also from the HS image, and can reduce the distortion caused by the insufficient spatial information.

The remainder of this paper is organized as follows. Section 2 describes the weighted structure tensor matrix and homomorphic filtering. In Section 3, the proposed homomorphic filtering and weighted tensor matrix-based hyperspectral pansharpening algorithm is presented. Experimental results and discussion are provided in Section 4, and conclusions are drawn in Section 5.

For an image $I$, the structure tensor matrix $M$ can be decomposed as:
where ${I}_{x}=\partial I/\partial x$ and ${I}_{y}=\partial I/\partial y$ are the horizontal and vertical partial derivatives of the image, $\nabla I={[\begin{array}{cc}{I}_{x}& {I}_{y}\end{array}]}^{\mathrm{T}}$, ${()}^{\mathrm{T}}$ is the transpose operation, ${v}_{1}$, ${v}_{2}$, ${e}_{1}$ and ${e}_{2}$ are the two eigenvectors and the corresponding eigenvalues, respectively. As shown in Equation (1), the tensor matrix $M$ which is a symmetric and semi-definite positive matrix has eigen-decomposition, and it has been exploited in some fields, such as texture synthesis [37], image regularization [38], denoising [39], and recognition systems [40]. The eigenvalues obtained by decomposing the tensor matrix are utilized to describe the structure information of the image. For the multiple images $[{I}_{1},{I}_{2},\dots ,{I}_{n}]$, the structural information is coupled from all images by employing the linear combination:
where ${M}_{w}$ is the weighted structure tensor matrix, ${M}_{l}$ is the structure tensor matrix of the $l\mathrm{th}$ image, $\partial {I}_{l}/\partial x$ and $\partial {I}_{l}/\partial y$ are the horizontal and vertical partial derivatives of the $l\mathrm{th}$ image, respectively. The weighted tensor matrix ${M}_{w}$ is also a symmetric and semi-definite positive matrix, and can be proceeded the eigen-decomposition.

$$M=\nabla I\xb7\nabla {I}^{\mathrm{T}}=\left[\begin{array}{cc}{I}_{x}^{2}& {I}_{x}{I}_{y}\\ {I}_{x}{I}_{y}& {I}_{y}^{2}\end{array}\right]=\left[\begin{array}{cc}{e}_{1}& {e}_{2}\end{array}\right]\left[\begin{array}{cc}{v}_{1}& 0\\ 0& {v}_{2}\end{array}\right]{\left[\begin{array}{cc}{e}_{1}& {e}_{2}\end{array}\right]}^{\mathrm{T}}$$

$${M}_{w}=\frac{1}{n}{\displaystyle {\sum}_{l=1}^{n}{M}_{l}}=\frac{1}{n}\left[\begin{array}{cc}{\displaystyle {\sum}_{l=1}^{n}{(\partial {I}_{l}/\partial x)}^{2}}& {\displaystyle {\sum}_{l=1}^{n}((\partial {I}_{l}/\partial x)(\partial {I}_{l}/\partial y))}\\ {\displaystyle {\sum}_{l=1}^{n}((\partial {I}_{l}/\partial x)(\partial {I}_{l}/\partial y))}& {\displaystyle {\sum}_{l=1}^{n}{(\partial {I}_{l}/\partial y)}^{2}}\end{array}\right]$$

Homomorphic filtering which is a type of frequency domain filtering can compress the image brightness range and enhance the image contrast. Homomorphic filtering has been applied to some image processing problems [41,42,43] and is based on an image imaging model:
where $f$ represents an image, ${f}_{\mathrm{H}}$ represents the high frequency reflectance component, and ${f}_{\mathrm{L}}$ represents the low frequency illumination component. Homomorphic filtering aims to reduce the low frequency component of an image. Logarithmic transformation is utilized to separate the two components:
After applying the Fourier transform:
where $F$, ${F}_{\mathrm{H}}$ and ${F}_{\mathrm{L}}$ denote the Fourier transform of $\mathrm{ln}(f)$, $\mathrm{ln}({f}_{\mathrm{H}})$ and $\mathrm{ln}({f}_{\mathrm{L}})$, respectively. Then, the high-pass filter $H$ is applied to Equation (5) as
where $S$ is the filtered result. The final image is obtained by the inverse Fourier transform and the exponential operation:
where ${f}_{hf}$ denotes the homomorphic filtered image, $s$ denotes the inverse Fourier transform of $S$, and ${\Im}^{-1}$ denotes the inverse Fourier transform.

$$f={f}_{\mathrm{H}}\xb7{f}_{\mathrm{L}}$$

$$\mathrm{ln}(f)=\mathrm{ln}({f}_{\mathrm{H}})+\mathrm{ln}({f}_{\mathrm{L}})$$

$$F={F}_{\mathrm{H}}+{F}_{\mathrm{L}}$$

$$S=F\xb7H={F}_{\mathrm{H}}\xb7H+{F}_{\mathrm{L}}\xb7H$$

$${f}_{hf}=\mathrm{exp}(s)=\mathrm{exp}({\Im}^{-1}(S))$$

Figure 2 shows the schematic of the proposed homomorphic filtering and weighted tensor (HFWT) matrix-based hyperspectral pansharpening algorithm. Let the LRHS image be represented by ${X}_{\mathrm{LR}}^{\mathrm{HS}}\in {\mathrm{R}}^{m\times n\times B}$, and the HRPAN image be denoted by ${X}^{\mathrm{PAN}}\in {\mathrm{R}}^{M\times N\times 1}$, where $m\times n$ and $M\times N$ are the size of the LRHS and the HRPAN images, respectively, and $B$ is the number of the LRHS image bands. The fused HRHS image is represented by ${X}_{\mathrm{HR}}^{\mathrm{HS}}\in {\mathrm{R}}^{M\times N\times B}$.

The open-closing operation which belongs to the mathematical morphology operation is an effective denoising processing operation [44,45]. The effects of denoising using the open operation and closed operation alone are usually not very good, since they may cause amplitude deflection. By contrast, the open-closing operation has the better denoising effect. The open operation is first applied on the image, and the selected structure element is larger than the noise size to remove the background noise. Then, the closed operation is utilized to remove the noise of the image obtained in the previous step. The open-closing denoising operation is suitable for the images which have less small details. Since the LRHS image has the low spatial resolution, the fine spatial details are few. The open-closing operation is applicable to removing noise with high interference in the LRHS image. The open-closing morphological operation is applied as:
for $k=1,2,\dots ,B$, where ${X}_{\mathrm{LR}}^{\mathrm{RNH}}$ denotes the denoised LRHS image, ${({X}_{\mathrm{LR}}^{\mathrm{HS}})}^{k}$ and ${({X}_{\mathrm{LR}}^{\mathrm{RNH}})}^{k}$ denote the $k\mathrm{th}$ band of the LRHS image and the denoised LRHS image, respectively, and ${S}_{1}$ and ${S}_{2}$ are the structure elements. Here, $\circ $ represents the opening operation which first uses the erosion operation and then the dilation operation, and $\u2022$ denotes the closing operation which does in reverse. The erosion and dilation operations obtain the local minimum and maximum of the image, respectively. Equation (8) can be expressed in detail as:
for $k=1,2,\dots ,B$, where $\mathsf{\Theta}$ and $\oplus $ denote the erosion and dilation operations, respectively.

$${({X}_{\mathrm{LR}}^{\mathrm{RNH}})}^{k}=({({X}_{\mathrm{LR}}^{\mathrm{HS}})}^{k}\circ {S}_{1})\u2022{S}_{2}$$

$${({X}_{\mathrm{LR}}^{\mathrm{OHS}})}^{k}={({X}_{\mathrm{LR}}^{\mathrm{HS}})}^{k}\circ {S}_{1}=({({X}_{\mathrm{LR}}^{\mathrm{HS}})}^{k}\mathsf{\Theta}{S}_{1})\oplus {S}_{1}=\underset{j\in {S}_{1}}{\mathrm{max}}\{[\underset{j\in {S}_{1}}{\mathrm{min}}({({X}_{\mathrm{LR}}^{\mathrm{HS}})}^{k}(i+j)-{S}_{1}(j))]+{S}_{1}(j)\}$$

$$({({X}_{\mathrm{LR}}^{\mathrm{HS}})}^{k}\circ {S}_{1}))\u2022{S}_{2}={({X}_{\mathrm{LR}}^{\mathrm{OHS}})}^{k}\u2022{S}_{2}=\underset{j\in {S}_{2}}{\mathrm{min}}\{[\underset{j\in {S}_{2}}{\mathrm{max}}({({X}_{\mathrm{LR}}^{\mathrm{OHS}})}^{k}(i+j)-{S}_{2}(j))]+{S}_{2}(j)\}$$

Homomorphic filtering is a filtering method that transforms the nonlinear problem into a linear problem. It transforms the nonlinear multiplicative mixed problem into an additive model by the logarithmic transformation, and then uses linear filtering to process it. The homomorphic filtering suppresses the low frequency illumination component and enhances the high frequency reflectance component. For an HS image, the high frequency component of each band is considered as the spatial component for each band. To obtain the spatial information for each band, we apply homomorphic filtering to each band of the denoised LRHS image. Through the use of homomorphic filtering, the low frequency component of each band of the denoised LRHS image is suppressed, and the high frequency component is extracted. Therefore, in this research, homomorphic filtering is applied to each band of the denoised LRHS image to extract the spatial component of each band. The homomorphic filtering processing is based on the following image imaging model:
for $k=1,2,\dots ,B$, where ${X}_{\mathrm{LR}\_\mathrm{H}}^{\mathrm{RNH}}$ represents the high frequency component of the denoised LRHS image, ${X}_{\mathrm{LR}\_\mathrm{L}}^{\mathrm{RNH}}$ represents the low frequency component, and ${({X}_{\mathrm{LR}\_\mathrm{H}}^{\mathrm{RNH}})}^{k}$ and ${({X}_{\mathrm{LR}\_\mathrm{L}}^{\mathrm{RNH}})}^{k}$ represent the $k\mathrm{th}$ band of ${X}_{\mathrm{LR}\_\mathrm{H}}^{\mathrm{RNH}}$ and ${X}_{\mathrm{LR}\_\mathrm{L}}^{\mathrm{RNH}}$, respectively. Based on Equation (4)-(6), Logarithmic transformation, Fourier transform, and high-pass filtering operations are applied to Equation (11):
for $k=1,2,\dots ,B$, where ${S}_{\mathrm{LR}}$ is the high-pass filtered imageh, ${({S}_{\mathrm{LR}})}^{k}$ is the $k\mathrm{th}$ band of ${S}_{\mathrm{LR}}$, $\Im $ represents Fourier transform, and $H$ is the high-pass filter, defined as:
where ${D}_{0}$ is the cut-off frequency, $D$ is the distance between $(x,y)$ and the center, ${\beta}_{\mathrm{H}}$ and ${\beta}_{\mathrm{L}}$ are the high and low frequency gains. Figure 3 shows the 3-D mesh of the high-pass filter. Since homomorphic filtering aims to reduce the low frequency component and extract the high frequency component, ${\beta}_{\mathrm{H}}$ is greater than 1 and ${\beta}_{\mathrm{L}}$ is smaller than 1. By adjusting the value of the cut-off frequency ${D}_{0}$, the sharpness of the transition between ${\beta}_{\mathrm{L}}$ and ${\beta}_{\mathrm{H}}$ can be controlled. In practice, the values of these parameters are generally determined empirically. In this paper, empirically, ${\beta}_{\mathrm{H}}$, ${\beta}_{\mathrm{L}}$, and ${D}_{0}$ are set to 2, 0.25, 40, respectively. ${S}_{\mathrm{LR}}$ is the high-pass filtered image in which the low frequency component has been weakened. Then, the spatial component of each band is obtained by applying the inverse Fourier transform and the exponential operation to ${S}_{\mathrm{LR}}$.
for $k=1,2,\dots ,B$, where ${X}_{\mathrm{LR}}^{\mathrm{I}}$ denotes the spatial component of each band of the denoised LRHS image, ${({X}_{\mathrm{LR}}^{\mathrm{I}})}^{k}$ denotes the $k\mathrm{th}$ band of ${X}_{\mathrm{LR}}^{\mathrm{I}}$, and ${\Im}^{-1}$ denotes the inverse Fourier transform.

$${({X}_{\mathrm{LR}}^{\mathrm{RNH}})}^{k}={({X}_{\mathrm{LR}\_\mathrm{H}}^{\mathrm{RNH}})}^{k}\xb7{({X}_{\mathrm{LR}\_\mathrm{L}}^{\mathrm{RNH}})}^{k}$$

$${({S}_{\mathrm{LR}})}^{k}=[\Im (\mathrm{ln}{({X}_{\mathrm{LR}\_\mathrm{H}}^{\mathrm{RNH}})}^{k})]\xb7H+\Im [\mathrm{ln}{({X}_{\mathrm{LR}\_\mathrm{L}}^{\mathrm{RNH}})}^{k}]\xb7H$$

$$H(x,y)=({\beta}_{\mathrm{H}}-{\beta}_{\mathrm{L}})[1-\mathrm{exp}(-({D}^{2}(x,y)/{D}_{0}^{2})]+{\beta}_{\mathrm{L}}$$

$${({X}_{\mathrm{LR}}^{\mathrm{I}})}^{k}=\mathrm{exp}[{\Im}^{-1}({({S}_{\mathrm{LR}})}^{k})]$$

After introducing homomorphic filtering to obtain the spatial component of each band of the denoised LRHS image, a weighted root mean squared error (RMSE)-based method is presented to extract the spatial intensity information of the HS image. Let ${I}_{\mathrm{LR}}={\displaystyle {\sum}_{k=1}^{B}{\lambda}_{k}}{({X}_{\mathrm{LR}}^{\mathrm{I}})}^{k}$ denote the total spatial information of the LRHS image, where $[{\lambda}_{1},{\lambda}_{2},\dots ,{\lambda}_{B}]$ is the weighted vector. To determine the values of the weighted vector, we utilize the RMSE index to measure the deviation of two images. A smaller value of RMSE indicates a better result, and the optimal value is 0. In the RMSE-based method, the spatial information of the PAN image is considered, and the RMSE value between the total spatial information ${I}_{\mathrm{LR}}$ and the PAN image ${X}^{\mathrm{PAN}}$ is calculated. The smallest value of RMSE is computed to obtain the optimal values of the weights $[{\lambda}_{1},{\lambda}_{2},\dots ,{\lambda}_{B}]$:
where $T=m\times n$ represents the total pixel number of one band of the LRHS image, $\downarrow $ represents down-sampling operation, $\downarrow {X}^{\mathrm{PAN}}$ represents that the PAN image is down-sampled to the size of the LRHS image, ${(\downarrow {X}^{\mathrm{PAN}})}_{i}$ and ${({\displaystyle {\sum}_{k=1}^{B}{\lambda}_{k}}{({X}_{\mathrm{LR}}^{\mathrm{I}})}^{k})}_{i}$ represent the values of the $i\mathrm{th}$ pixel in $\downarrow {X}^{\mathrm{PAN}}$ and ${\sum}_{k=1}^{B}{\lambda}_{k}}{({X}_{\mathrm{LR}}^{\mathrm{I}})}^{k$, respectively. The laplacian pyramid networks (LapSRN) [46] super-resolution method can effectively improve the spatial resolution of an image, and has the advantages of parameter sharing, local skip connections, and multi-scale training. So it is adopted to super-resolve the spatial information of the LRHS image ${I}_{\mathrm{LR}}$ for an ${I}_{\mathrm{HR}}$ with super-resolution spatial information.

$$\mathrm{min}(\sqrt{\frac{1}{T}{\displaystyle \sum _{i=1}^{T}{[{({\displaystyle {\sum}_{k=1}^{B}{\lambda}_{k}}{({X}_{\mathrm{LR}}^{\mathrm{I}})}^{k})}_{i}-{(\downarrow {X}^{\mathrm{PAN}})}_{i}]}^{2}}}$$

To make the spatial information of the PAN image clearer, the Laplacian of Gaussian (LOG) [47] image enhancement algorithm is applied to the PAN image, which uses a Gaussian filter to reduce noise followed by a Laplace operator for enhancement. Let ${I}_{\mathrm{s}}^{\mathrm{PAN}}$ represent the enhanced PAN image.

The HS and PAN images contain the different and complementary information for a scene. To acquire the total spatial information, we should consider simultaneously the spatial structure details of these two images, and then propose an optimized weighted tensor matrix-based method. ${I}_{\mathrm{HR}}$ and ${I}_{\mathrm{s}}^{\mathrm{PAN}}$ include the spatial structure information of the HS and PAN images, respectively. Based on Equation (2), for the multiple images $[{I}_{\mathrm{HR}},{I}_{\mathrm{s}}^{\mathrm{PAN}}]$, the weighted structure tensor matrix at pixel $p$ is given by:
where ${M}_{w}^{\mathrm{HP}}$ denotes the obtained weighted tensor matrix, ${M}_{w,p}^{\mathrm{HP}}$ denotes ${M}_{w}^{\mathrm{HP}}$ at pixel $p$, $\partial ({I}_{\mathrm{HR},p})/\partial x$, $\partial ({I}_{\mathrm{HR},p})/\partial y$, $\partial ({I}_{\mathrm{s},p}^{\mathrm{PAN}})/\partial x$, and $\partial ({I}_{\mathrm{s},p}^{\mathrm{PAN}})/\partial y$ are the $x$ and $y$ partial derivatives of ${I}_{\mathrm{HR}}$ and ${I}_{\mathrm{s}}^{\mathrm{PAN}}$ at pixel $p$, respectively. The weighted tensor matrix ${M}_{w,p}^{\mathrm{HP}}$ is semi-definite, and it can be decomposed as:
where ${()}^{\mathrm{T}}$ is the transpose operation, ${v}_{w1}$ and ${v}_{w2}$ are the two eigenvalues, ${v}_{w1,p}$ and ${v}_{w2,p}$ are the two eigenvalues at pixel $p$, ${e}_{w1,p}={[\begin{array}{cc}{e}_{w11,p}& {e}_{w12,p}\end{array}]}^{\mathrm{T}}$ and ${e}_{w2,p}={[\begin{array}{cc}{e}_{w21,p}& {e}_{w22,p}\end{array}]}^{\mathrm{T}}$ are the eigenvectors corresponding to the two eigenvalues at pixel $p$, respectively.

$${M}_{w,p}^{\mathrm{HP}}=\frac{1}{2}\left[\begin{array}{cc}{(\partial ({I}_{\mathrm{HR},p})/\partial x)}^{2}+{(\partial ({I}_{\mathrm{s},p}^{\mathrm{PAN}})/\partial x)}^{2}& (\partial ({I}_{\mathrm{HR},p})/\partial x)(\partial ({I}_{\mathrm{HR},p})/\partial y)+(\partial ({I}_{\mathrm{s},p}^{\mathrm{PAN}})/\partial x)(\partial ({I}_{\mathrm{s},p}^{\mathrm{PAN}})/\partial y)\\ (\partial ({I}_{\mathrm{HR},p})/\partial x)(\partial ({I}_{\mathrm{HR},p})/\partial y)+(\partial ({I}_{\mathrm{s},p}^{\mathrm{PAN}})/\partial x)(\partial ({I}_{\mathrm{s},p}^{\mathrm{PAN}})/\partial y)& {(\partial ({I}_{\mathrm{HR},p})/\partial y)}^{2}+{(\partial ({I}_{\mathrm{s},p}^{\mathrm{PAN}})/\partial y)}^{2}\end{array}\right]$$

$$\begin{array}{cc}\hfill {M}_{w,p}^{\mathrm{HP}}& =\left[\begin{array}{cc}{e}_{w11,p}& {e}_{w21,p}\\ {e}_{w12,p}& {e}_{w22,p}\end{array}\right]\left[\begin{array}{cc}{v}_{w1,p}& 0\\ 0& {v}_{w2,p}\end{array}\right]{\left[\begin{array}{cc}{e}_{w11,p}& {e}_{w21,p}\\ {e}_{w12,p}& {e}_{w22,p}\end{array}\right]}^{\mathrm{T}}\hfill \\ & ={v}_{w1,p}{e}_{w1,p}{({e}_{w1,p})}^{\mathrm{T}}+{v}_{w2,p}{e}_{w2,p}{({e}_{w2,p})}^{\mathrm{T}}\hfill \end{array}$$

The two eigenvalues generally have a larger value and a smaller value. We assume that ${v}_{w1}$ is the larger eigenvalue. When ${v}_{w1}\approx {v}_{w2}\approx 0$, ${v}_{w1}>{v}_{w2}\approx 0$, and ${v}_{w1}>{v}_{w2}>0$, the structure region of this pixel are flat area, edge area, and corner, respectively. We test on some images to study the eigenvalues of weighted tensor matrix. Figure 4 shows the two eigenvalues at each pixel of the weighted tensor matrix for the Salinas scene data. It can be seen that for many pixels, the smaller eigenvalues shown in Figure 4d are approximately ${10}^{-5}$, and are very small. By experimenting on other numerous images, we have also discovered that the smaller eigenvalues are mostly very small. Thus, the approximation of ${M}_{w,p}^{\mathrm{HP}}$ is expressed as:
where ${\tilde{M}}_{w,p}^{\mathrm{HP}}$ is the approximation of ${M}_{w,p}^{\mathrm{HP}}$. Based on Equation (1), the structure tensor matrix satisfies that $M=\nabla I\xb7\nabla {I}^{\mathrm{T}}$. The weighted gradient ${G}_{w}$ at pixel $p$ satisfies that ${\tilde{M}}_{w,p}^{\mathrm{HP}}={G}_{w,p}\xb7{({G}_{w,p})}^{\mathrm{T}}$$={v}_{w1,p}{e}_{w1,p}{({e}_{w1,p})}^{\mathrm{T}}$. So, ${G}_{w,p}$ is deduced as ${G}_{w,p}=\sqrt{{v}_{w1,p}}\xb7{e}_{w1,p}$. Since the direction of the eigenvector corresponding to ${v}_{w1}$ is not unique, the direction of the weighted gradient ${G}_{w,p}$ is also not unique. We specify the direction of the weighted gradient ${G}_{w,p}$ as the gradients average of the individual multiple source images $[{I}_{\mathrm{HR}},{I}_{\mathrm{s}}^{\mathrm{PAN}}]$:
where $\nabla {I}_{\mathrm{HR},p}=[\begin{array}{cc}\partial ({I}_{\mathrm{HR},p})/\partial x& \partial ({I}_{\mathrm{HR},p})/\partial y\end{array}]$ and $\nabla {I}_{\mathrm{s},p}^{\mathrm{PAN}}=[\begin{array}{cc}\partial ({I}_{\mathrm{s},p}^{\mathrm{PAN}})/\partial x& \partial ({I}_{\mathrm{s},p}^{\mathrm{PAN}})/\partial y\end{array}]$, $\langle \cdot ,\cdot \rangle $ represents the inner product of two vectors, and $\mathrm{sign}(\cdot )$ represents the sign function. Once the weighted gradient ${G}_{w,p}$ is acquired from the multiple images $[{I}_{\mathrm{HR}},{I}_{\mathrm{s}}^{\mathrm{PAN}}]$, an optimization model is proposed to obtain the total spatial information ${I}_{\mathrm{T}}^{\mathrm{HP}}$ as:
where $\nabla {I}_{\mathrm{T}}^{\mathrm{HP}}=[\begin{array}{cc}\partial {I}_{\mathrm{T}}^{\mathrm{HP}}/\partial x& \partial {I}_{\mathrm{T}}^{\mathrm{HP}}/\partial y\end{array}]$, $\partial {I}_{\mathrm{T}}^{\mathrm{HP}}/\partial x$, and $\partial {I}_{\mathrm{T}}^{\mathrm{HP}}/\partial y$ denote the $x$ and $y$ partial derivatives of ${I}_{\mathrm{T}}^{\mathrm{HP}}$. Equation (20) is an unconstrained optimization problem, and we solve it by the conjugate gradient method. Equation (20) can effectively ensure that the total spatial information ${I}_{\mathrm{T}}^{\mathrm{HP}}$ contains the spatial structure details of both the HS and PAN images.

$${\tilde{M}}_{w,p}^{\mathrm{HP}}=\left[\begin{array}{cc}{e}_{w11,p}& {e}_{w21,p}\\ {e}_{w12,p}& {e}_{w22,p}\end{array}\right]\left[\begin{array}{cc}{v}_{w1,p}& 0\\ 0& 0\end{array}\right]{\left[\begin{array}{cc}{e}_{w11,p}& {e}_{w21,p}\\ {e}_{w12,p}& {e}_{w22,p}\end{array}\right]}^{\mathrm{T}}={v}_{w1,p}{e}_{w1,p}{({e}_{w1,p})}^{\mathrm{T}}$$

$${G}_{w,p}=\sqrt{{v}_{w1,p}}\xb7{e}_{w1,p}\xb7\mathrm{sign}\langle {e}_{w1,p},\frac{1}{2}(\nabla {I}_{\mathrm{HR},p}+\nabla {I}_{\mathrm{s},p}^{\mathrm{PAN}})\rangle $$

$$\mathrm{min}{\Vert (\nabla {I}_{\mathrm{T}}^{\mathrm{HP}})-{G}_{w}\Vert}^{2}$$

The LRHS image ${X}_{\mathrm{LR}}^{\mathrm{HS}}$ is interpolated to the scale of the HRPAN image. By constructing a suitable gain matrix $R$, the total spatial information ${I}_{\mathrm{T}}^{\mathrm{HP}}$ is injected into the interpolated HS image to generate the fused HRHS image ${X}_{\mathrm{HR}}^{\mathrm{HS}}$. For the gain matrix $R$, it is beneficial to preserve the ratio of each HS pair band unchanged to reduce the spectral distortion. Thus, $R$ should satisfy that ${R}^{k}\propto \frac{{({X}_{\mathrm{IN}}^{\mathrm{HS}})}^{k}}{(1/B){\displaystyle {\sum}_{k=1}^{B}{({X}_{\mathrm{IN}}^{\mathrm{HS}})}^{k}}}$, where ${X}_{\mathrm{IN}}^{\mathrm{HS}}$ is the interpolated HS image, ${({X}_{\mathrm{HR}}^{\mathrm{HS}})}^{k}$ and ${R}^{k}$ are the $k\mathrm{th}$ band of ${X}_{\mathrm{IN}}^{\mathrm{HS}}$ and $R$, respectively. Then, a tradeoff parameter $\epsilon $ is defined to regulate the amount of the injected details to reduce the spatial distortion. This process can be expressed as:
for $k=1,2,\dots ,B$, where ${X}_{\mathrm{HR}}^{\mathrm{HS}}$ is the fused HRHS image, and ${({X}_{\mathrm{HR}}^{\mathrm{HS}})}^{k}$ is the $k\mathrm{th}$ band of ${X}_{\mathrm{HR}}^{\mathrm{HS}}$.

$${({X}_{\mathrm{HR}}^{\mathrm{HS}})}^{k}={({X}_{\mathrm{IN}}^{\mathrm{HS}})}^{k}+{R}^{k}\xb7{I}_{\mathrm{T}}^{\mathrm{HP}}={({X}_{\mathrm{IN}}^{\mathrm{HS}})}^{k}+\epsilon \xb7\frac{{({X}_{\mathrm{IN}}^{\mathrm{HS}})}^{k}}{(1/B){\displaystyle {\sum}_{k=1}^{B}{({X}_{\mathrm{IN}}^{\mathrm{HS}})}^{k}}}\xb7{I}_{\mathrm{T}}^{\mathrm{HP}}$$

In order to evaluate the effectiveness of the proposed HFWT hyperspectral pansharpening method (named as HFWT), experiments were performed on two simulated hyperspectral datasets which were a Washington DC and a Salinas scene, and one real hyperspectral dataset, the Hyperion dataset. The Salinas scene hyperspectral dataset was collected by Airborne Visible Infrared Imaging Spectrometer (AVIRIS) [48], and the Washington DC dataset was acquired by the Spectral Information Technology Application Center of Virginia. The used real dataset is provided by the EO-1 spacecraft. The EO-1 spacecraft has a Hyperion instrument which provides the real LRSH images and an Advanced Land Imager (ALI) instrument acquires the HRPAN images [48]. Table 1 lists the characteristics of each dataset.

The proposed HFWT method is compared with several state-of-the-art hyperspectral pansharpening methods: Gram-Schmidt (GS) [23], guided filter principal component analysis (GFPCA) [33], coupled nonnegative matrix factorization (CNMF) [18], Bayesian sparsity promoted Gaussian prior (Bayesian) [13], and HySure [14]. Four typical quantitative evaluation indexes are adopted: cross correlation (CC) [49], spectral angle mapper (SAM) [50], root mean squared error (RMSE), and erreur relative globale adimensionnelle desynthèse (ERGAS) [51]. The CC and SAM measure the spectral and spatial distortion, respectively. The larger value of CC and the smaller value of SAM indicate the better fusion result. The RMSE and ERGAS are the global indexes that measure both the spatial and spectral performance, and their value ranges are all (0,1), with 0 being the optimal value.

In order to perform the objective fusion evaluation of the simulated hyperspectral datasets, the available HS image is used as the reference HS image. The simulated LRHS image and HRPAN image are generated according to Wald’s protocol [52,53]. The reference HS image is blurred and down sampled 4 times to obtain the simulated HS image. The simulated PAN image is obtained by averaging the visible light band of the reference HS image.

For the real datasets, the real LRHS and HRPAN images are available, and the reference high resolution HS image is not available. In order to test the objective quality of the real hyperspectral images, the real LRHS image is served as the reference image. The real LRHS and HRPAN images available are degraded, and the two degraded images are fused to obtain a fused image. This fused image is compared to the real LRHS image to evaluate the objective quality.

In the proposed HFWT method, we define a tradeoff parameter $\epsilon $ which regulates the amount of the injected details and ensures spatial performance. In practice, the optimal value is determined based on experience. By adjusting different values of $\epsilon $, the optimal value can be determined by the fusion result. In this paper, by experience, the values of tradeoff parameter $\epsilon $ are set as 0.25, 0.05 and 0.2 for the Washington DC, Salinas scene and Hyperion dataset, respectively.

To verify the effectiveness of the open-closing HS image denoising operation, the proposed HFWT method was conducted on the Washington DC dataset with different denoising processing. The compared denoising algorithms contain average filtering, Gaussian filtering, open operation, closed operation, and open-closing operation. Table 2 shows the fusion performance of different image denoising processing. As outlined in Table 2, the HFWT method without HS image denoising had the worst fusion results. The proposed methods with each HS image denoising preprocessing have the better fusion results compared with the proposed method without HS image denoising, which demonstrates that the HS image denoising preprocessing is significative and effective. By contrast, the proposed HFWT method with the open-closing operation achieves the best fusion performance, and it demonstrates that the open-closing operation is an effective HS image denoising preprocessing.

Figure 5 shows the fusion experimental results for the Washington DC dataset, where Figure 5(a1) shows the reference HS image, and the subjective fused HS images of each method are displayed in Figure 5(b1–g1). Moreover, Figure 5(a2) shows the enlarged subareas of the reference HS image, and the two enlarged subareas of each fused image are shown in Figure 5(b2–g2). The reference error image is shown in Figure 5(a3), and the error images between each fused HS image and the reference HS image are reported in Figure 5(b3–g3). Except for the first column, each column in Figure 5 shows the experimental results corresponding to each method. By visually comparing the fused HS images with the reference HS image, the fused result of the GS method suffers from serious spectral distortion. For example, the two enlarged subareas of the GS method are distorted seriously. The GFPCA approach generates fuzzy spatial details in some regions, such as the enlarged subareas shown in Figure 5(c2). This is because the spatial information of the GFPCA approach is injected insufficiently. As depicted in Figure 5(d1,d2), the spatial information of the fused image is well enhanced using the CNMF method, but some slight spectral distortion is appeared in the roof of the buildings. A closer inspection revealed that the HySure method seems to generate some distortion in the circular building in the upper left corner. By contrast, the fused HS images obtained by the Bayesian and HFWT methods achieve superior performance in terms of both spectral and spatial aspects. In order to further compare the performance of each fusion method, the third row of Figure 5 shows the error images of different methods. The error image is the difference (absolute value) of pixel values between the fused HS images and the reference HS image. We can see that, the GS, GFPCA and CNMF methods have larger differences, the HySure and Bayesian approaches generated relatively smaller differences, and the proposed HFWT approach shows the smallest differences in most areas, which demonstrates the excellent fusion capacity of the proposed method.

Similar to the previous experiments, for the Salinas scene dataset, the fused results are shown in Figure 6. Figure 6(a1–a3) show the reference HS image, the enlarged subareas of the reference HS image, and the reference SAM image, respectively. Figure 6(b1–g1) in the first row show the subjective pansharpened results of each algorithm, and Figure 6(b2–g2) in the second row display each enlarged subarea. The SAM images of each approach are shown in Figure 6(b3–g3). The reference SAM image of enlarged subarea is shown in Figure 6(a4), and the SAM images of enlarged subarea obtained by each method are shown in Figure 6(b4–g4). The spectral distortion caused by the GS method is very obvious, and the degree of spatial enhancement is also not acceptable for the GS method, as depicted on Figure 6(b1,b2). Compared with the GS approach, the GFPCA method performs better in terms of the spectral quality. However, the fused HS image obtained by the GFPCA approach shows an indistinct area in the left region of Figure 6(c1). Despite having a preeminent spatial quality, the CNMF method generates significant spectral distortion in the triangle region in the lower half of Figure 6(d1). From the visual analysis, the HySure, Bayesian and HFWT methods effectively improve spatial performance while maintaining spectral information, and the HFWT method shows better spectral quality compared with the HySure and Bayesian methods in some regions, such as the upper area of the enlarged subarea. The SAM images and the SAM images of the enlarged subareas of different approaches are shown in the third and fourth rows of Figure 6, to further verify the fusion performance of the proposed method. It can be seen that the proposed HFWT method yields the lowest SAM values for most regions. These results demonstrate that the proposed HFWT algorithm performs well in both the spatial and spectral aspects.

In addition to visual inspection, the performance of each algorithm for the Washington DC and Salinas scene datasets is analyzed quantitatively in Table 3, where the best results for each quantitative index are marked in bold. As can be seen from Table 3, the objective quantitative results are roughly consistent with the subjective qualitative effects. Same as the subjective results, the GS and GFPCA algorithms produce worse objective performance compared with other algorithms. The HySure approach obtains the best RMSE value for the Washington DC dataset and the optimal ERGAS value for the Salinas scene dataset. Most of the quality indexes generated by applying the proposed HFWT method are the best, in which the SAM, CC and ERGAS values are the best for Washington DC, and the RMSE, SAM and CC indexes are ranked first for the Salinas scene.

Figure 7 shows the pansharpened images of each method for the Hyperion dataset to confirm the fusion performance of the proposed HFWT method in the real dataset. Figure 7a–c show the real HS, real PAN, and interpolated HS images, respectively. The GS method shown in Figure 7d generates obvious spectral distortion, especially in the wharf area. In spite of good spatial improvement, the spatial details shown in the fused images obtained by using the GS and HySure approaches are too sharp. Spectral quality of the GFPCA and Bayesian methods seems be acceptable, but the GFPCA and Bayesian methods perform poorly from the spatial aspect. By contrast, the subjective effects of the CNMF and HFWT approaches are the best, and the HFWT method yields better spatial capacity compared to the CNMF method. The objective quality evaluation for the Hyperion dataset are presented in Table 4. As reported in Table 4, the HFWT method provides the best quantitative evaluation results in terms of the RMSE, SAM, CC, and ERGAS indices, which indicates that the HFWT method successfully maintains the spectral information of the original LRHS image and improves the spatial resolution.

The proposed HFWT algorithm contains simple sequential statements, several loop statements without nesting, and a two-layer loop statements with nesting (In the proposed HFWT algorithm, the program statement of Equation (19) which are applied on each pixel is a two-layer loop statements). The simple sequential statement is naturally $O(1)$ time, and the loop statement without nesting is $O(n)$ time. The two-layer loop statement is $O({n}^{2})$. According to the summation rule of algorithm complexity, the total algorithm complexity is $O({n}^{2}+n+1)=O({n}^{2})$. The proposed HFWT algorithm which is $O({n}^{2})$ time belongs to the polynomial time, and is considered a fast algorithm. The computing time (in seconds) of each method for three datasets is shown in Table 5. The experiments in this paper were all performed using MATLAB R2015b, and tested on a PC with an Intel Core i5-7300HQ CPU @ 2.50 GHz and 8 GB of memory. The GS and GFPCA methods are very efficient, but the fusion performance of the GS and GFPCA methods is unsatisfactory. The proposed HFWT method is faster than the CNMF algorithm, and takes much less computing time than the HySure and Bayesian algorithms. The time cost of the proposed HFWT is acceptable.

This paper presents a novel hyperspectral pansharpening method based on the merger of the homomorphic filtering and weighted tensor matrix. The proposed HFWT algorithm introduces the open-closing morphological operation and homomorphic filtering to remove noise and extract spatial information of each band of an HhS image, respectively. Moreover, we propose a weighted RMSE-based method to obtain the total spatial information of the HS image. In order to generate the adequate spatial information from both the HS and the corresponding PAN images, an optimized weighted tensor matrix based method is proposed. Specifically, the weighted tensor matrix, eigenvalues and eigenvectors are deduced and analyzed to obtain the weighted gradient, and an optimization model is presented to acquire the integrated spatial information. Compared with the state-of-the-art methods, experiments performed on the Washington DC, Salinas scene and Hyperion datasets demonstrate the proposed method performs superiorly in terms of both subjective and objective assessment.

Methodology, J.Q., Y.L. and Q.D.; Software, J.Q. and W.D.; Supervision, Y.L. and Q.D.; Writing—original draft, J.Q.; Writing—review and editing, Y.L., Q.D., W.D. and B.X.

This work was supported in part by the National Natural Science Foundation of China (nos. 61571345, 91538101, 61501346, 61,502,367 and 61701360) and the 111 project (B08038). It was also partially supported by the Supported by Yangtze River Scholar Bonus Schemes of China (No. CJT160102), Ten Thousand Talent Program, the Natural Science Basic Research Plan in Shaanxi Province of China (No. 2016JQ6023), the China Scholarship Council program (201806960020), the Excellent doctoral thesis fund of Xidian University, the Innovation Fund of Xidian University, and Fundamental Research Funds for Central Universities (JB182001).

The authors would like to thank the editors and the anonymous reviewers for their insightful comments and suggestions which have greatly improved this paper.

The authors declare no conflict of interest.

- Adam, E.; Mutanga, O.; Rugege, D. Multispectral and hyperspectral remote sensing for identification and mapping of wetland vegetation: A review. Wetl. Ecol. and Manag.
**2010**, 18, 281–296. [Google Scholar] [CrossRef] - Zhou, Y.; Feng, L.Y.; Hou, C.P.; Kung, S.Y. Hyperspectral and multispectral image fusion based on local low rank and coupled spectral unmixing. IEEE Trans. Geosci. Remote Sens.
**2017**, 55, 5997–6009. [Google Scholar] [CrossRef] - Xu, Y.; Wu, Z.; Li, J.; Plaza, A.; Wei, Z. Anomaly detection in hyperspectral images based on low-rank and sparse representation. IEEE Trans. Geosci. Remote Sens.
**2016**, 54, 1990–2000. [Google Scholar] [CrossRef] - Zhang, M.; Li, W.; Du, Q. Diverse region-based CNN for hyperspectral image classification. IEEE Trans. Image Process.
**2018**, 27, 2623–2634. [Google Scholar] [CrossRef] - Xue, B.; Yu, C.; Wang, Y.; Song, M.; Li, S.; Wang, L.; Chen, H.; Chang, C.I. A subpixel target detection approach to hyperspectral image classification. IEEE Trans. Geosci. Remote Sens.
**2017**, 55, 5093–5114. [Google Scholar] [CrossRef] - Kang, X.; Duan, P.; Xiang, X.; Li, S.; Benediktsson, J.A. Detection and correction of mislabeled training samples for hyperspectral image classification. IEEE Trans. Geosci. Remote Sens.
**2018**, 56, 5673–5686. [Google Scholar] [CrossRef] - Lv, Z.Y.; Shi, W.Z.; Zhou, X.C.; Benediktsson, J.A. Semi-automatic system for land cover change detection using bi-temporal remote sensing images. Remote Sens.
**2017**, 9, 1112. [Google Scholar] [CrossRef] - Haboudane, D.; Miller, J.R.; Pattey, E.; Zarco-Tejada, P.; Ian, B.; Strachan, I. Hyperspectral vegetation indices and novel algorithms for predicting green LAI of crop canopies: Modeling and validation in the context of precision agriculture. Remote Sens. Environ.
**2004**, 90, 337–352. [Google Scholar] [CrossRef] - Chabrillat, S.; Pinet, P.C.; Ceuleneer, G.; Johnson, P.E.; Mustard, J.F. Ronda peridotite massif: Methodology for its geological mapping and lithological discrimination from airborne hyperspectral data. Int. J. Remote Sens.
**2000**, 21, 2363–2388. [Google Scholar] [CrossRef] - Bishop, C.A.; Liu, J.G.; Mason, P.J. Hyperspectral remote sensing for mineral exploration in Pulang, Yunnan Province, China. Int. J. Remote Sens.
**2011**, 32, 2409–2426. [Google Scholar] [CrossRef] - Ellis, R.J.; Scott, P.W. Evaluation of hyperspectral remote sensing as a means of environmental monitoring in the St. Austell China clay (kaolin) region, Cornwall, UK. Remote Sens. Environ.
**2004**, 93, 118–130. [Google Scholar] [CrossRef] - Wei, Q.; Dobigeon, N.; Tourneret, J.Y. Fast fusion of multiband images based on solving a Sylvester equation. IEEE Trans. Image Process.
**2015**, 24, 4109–4121. [Google Scholar] [CrossRef] - Wei, Q.; Dobigeon, N.; Tourneret, J.Y. Bayesian fusion of multiband images. IEEE J. Sel. Top. Signal Process.
**2015**, 9, 1117–1127. [Google Scholar] [CrossRef] - Simoes, M.; Dias, J.B.; Almeida, L.; Chanussot, J. A convex formulation for hyperspectral image superresolution via subspace-based regularization. IEEE Trans. Geosci. Remote Sens.
**2015**, 53, 3373–3388. [Google Scholar] [CrossRef] - Loncan, L.; Almeida, L.; Dias, J.B.; Briottet, X.; Chanussot, J.; Dobigeon, N.; Fabre, S.; Liao, W.Z.; Licciardi, G.A.; Simoes, M.; Tourneret, J.; Veganzones, M.A.; Vivone, G.; Wei, Q.; Yokoya, N. Hyperspectral pansharpening: A review. IEEE Geosci. Remote Sens. Mag.
**2015**, 3, 27–46. [Google Scholar] [CrossRef] - Yokoya, N.; Grohnfeldt, C.; Chanussot, J. Hyperspectral and multispectral data fusion: A comparative review of the recent literature. IEEE Geosci. Remote Sens. Mag.
**2017**, 5, 29–56. [Google Scholar] [CrossRef] - Hoyer, P.O. Non negative sparse coding. In Proceedings of the IEEE Workshop Neural Network Signal Processing, Martigny, Switzerland, 6 September 2002; pp. 557–565. [Google Scholar]
- Yokoya, N.; Yairi, T.; Iwasaki, A. Coupled nonnegative matrix factorization unmixing for hyper-spectral and multispectral data fusion. IEEE Trans. Geosci. Remote Sens.
**2012**, 50, 528–537. [Google Scholar] [CrossRef] - Carper, W.; Lillesand, T.M.; Kiefer, P.W. The use of Intensity-Hue-Saturation transformations for merging SPOT panchromatic and multispectral image data. Photogramm. Eng. Remote Sens.
**1990**, 56, 459–467. [Google Scholar] - Tu, T.M.; Su, S.C.; Shyu, H.C.; Huang, P.S. A new look at IHS-like image fusion methods. Inf. Fusion.
**2001**, 2, 117–186. [Google Scholar] [CrossRef] - Chavez, P.S.; Kwarteng, A.Y.A. Extracting spectral contrast in Landsat thematic mapper image data using selective principal component analysis. Photogramm. Eng. Remote Sens.
**1989**, 55, 339–348. [Google Scholar] - Shettigara, V. A generalized component substitution technique for spatial enhancement of multispectral images using a higher resolution data set. Photogramm. Eng. Remote Sens.
**1992**, 58, 561–567. [Google Scholar] - Laben, C.; Brower, B. Process for Enhacing the Spatial Resolution of Multispectral Imagery Using Pan-Sharpening. U.S. Patent 6,011,875, 4 January 2000. [Google Scholar]
- Aiazzi, B.; Baronti, S.; Selva, M. Improving component substitution pansharpening through multivariate regression of MS + Pan data. IEEE Trans. Geosci. Remote Sens.
**2007**, 45, 3230–3239. [Google Scholar] [CrossRef] - Thomas, C.; Ranchin, T.; Wald, L.; Chanussot, J. Synthesis of multispectral images to high spatial resolution: A critical review of fusion methods based on remote sensing physics. IEEE Trans. Geosci. Remote Sens.
**2008**, 46, 1301–1312. [Google Scholar] [CrossRef] - Liu, J.G. Smoothing filter based intensity modulation: A spectral preserve image fusion technique for improving spatial details. Int. J. Remote Sens.
**2000**, 21, 3461–3472. [Google Scholar] [CrossRef] - Burt, P.J.; Adelson, E.H. The Laplacian pyramid as a compact image code. IEEE Trans. Commun.
**1983**, 31, 532–540. [Google Scholar] [CrossRef] - Aiazzi, B.; Alparone, L.; Baronti, S.; Garzelli, A.; Selva, M. MTF-tailored multiscale fusion of high-resolution MS and pan imagery. Photogramm. Eng. Remote Sens.
**2006**, 72, 591–596. [Google Scholar] [CrossRef] - Vivone, G.; Restaino, R.; Mura, M.D.; Licciardi, G.; Chanussot, J. Contrast and error-based fusion schemes for multispectral image pansharpening. IEEE Geosci. Remote Sens. Lett.
**2014**, 11, 930–934. [Google Scholar] [CrossRef] - Aiazzi, B.; Alparone, L.; Baronti, S.; Garzelli, A.; Selva, M. 25 years of pansharpening: A critical review and new developments. In Signal Image Processing for Remote Sensing, 2nd ed.; Chen, C.H., Ed.; CRC Press: Boca Raton, FL, USA, 2011; Chapter 28; pp. 533–548. [Google Scholar]
- Park, H.; Choi, J.; Park, N.; Choi, S. Sharpening the VNIR and SWIR bands of Sentinel-2A imagery through modified selected and synthesized band schemes. Remote Sens.
**2017**, 9, 80. [Google Scholar] [CrossRef] - Vivone, G.; Addesso, P.; Restaino, R.; Dalla, M.; Chanussot, J. Pansharpening Based on Deconvolution for Multiband Filter Estimation. IEEE Trans. Geosci. Remote Sens.
**2019**, 57, 540–553. [Google Scholar] [CrossRef] - Liao, W.; Huang, X.; Coillie, F.; Gautama, S.; Pizurica, A.; Philips, W.; Liu, H.; Zhu, T.; Shimoni, M.; Moser, G.; et al. Processing of multiresolution thermal hyperspectral and digital color data: Outcome of the 2014 IEEE GRSS data fusion contest. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.
**2015**, 8, 2984–2996. [Google Scholar] [CrossRef] - Yang, J.; Zhao, Y.Q.; Chan, J.C.W. Hyperspectral and Multispectral Image Fusion via Deep Two-Branches Convolutional Neural Network. Remote Sens.
**2018**, 10, 800. [Google Scholar] [CrossRef] - Zhang, Y.J.; Liu, C.; Sun, M.W.; Ou, Y.J. Pan-Sharpening Using an Efficient Bidirectional Pyramid Network. IEEE Trans. Geosci. Remote Sens.
**2019**, 1–15. [Google Scholar] [CrossRef] - Palsson, F.; Sveinsson, J.R.; Ulfarsson, M.O. Multispectral and Hyperspectral Image Fusion Using a 3-D-Convolutional Neural Network. IEEE Geosci. Remote Sens. Lett.
**2017**, 14, 639–642. [Google Scholar] [CrossRef] - Akl, A.; Yaacoub, C.; Donias, M.; Costa, J.D.; Germain, C. Texture synthesis using the structure tensor. IEEE Trans. Image Process.
**2015**, 24, 4082–4095. [Google Scholar] [CrossRef] - Lefkimmiatis, S.; Osher, S. Nonlocal structure tensor functionals for image regularization. IEEE Trans. Computat. Imag.
**2015**, 1, 16–29. [Google Scholar] [CrossRef] - Wu, Z.; Wang, Q.; Jin, J.; Shen, Y. Structure tensor total variation-regularized weighted nuclear norm minimization for hyperspectral image mixed denoising. Signal Process.
**2017**, 131, 202–219. [Google Scholar] [CrossRef] - Tiwari, K.; Arya, D.K.; Badrinath, G.S.; Gupta, P. Designing palmprint based recognition system using local structure tensor and force field transformation for human identification. Neurocomputing.
**2013**, 116, 222–230. [Google Scholar] [CrossRef] - Fan, C.N.; Zhang, F.Y. Homomorphic filtering based illumination normalization method for face recognition. Pattern Recogn. Lett.
**2011**, 32, 1468–1479. [Google Scholar] [CrossRef] - Sreenivasan, K.R.; Havlicek, M.; Deshpande, G. Nonparametric hemodynamic deconvolution of fMRI using homomorphic filtering. IEEE Trans. Med. Imaging.
**2014**, 34, 1155–1163. [Google Scholar] [CrossRef] [PubMed] - Xiao, L.; Li, C.; Wu, Z.; Wang, T. An enhancement method for X-ray image via fuzzy noise removal and homomorphic filtering. Neurocomputing.
**2016**, 195, 56–64. [Google Scholar] [CrossRef] - Liu, Z.G.; Wang, J.L.; Liu, B. ECG Signal Denoising Based on Morphological Filtering. In Proceedings of the 2011 5th International Conference on Bioinformatics and Biomedical Engineering, Wuhan, China, 10–12 May 2011. 12072479. [Google Scholar]
- Chen, S.H.; Wang, E.Y. Electromagnetic Radiation Signals of Coal or Rock Denoising Based on Morphological Filter. Procedia Engineer.
**2011**, 26, 588–594. [Google Scholar] - Lai, W.; Huang, J.; Ahuja, N.; Yang, M. Deep laplacian pyramid networks for fast and accurate super-resolution. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Honolulu, HI, USA, 21–26 July 2017; pp. 624–632. [Google Scholar]
- Dong, W.Q.; Xiao, S.; Li, Y.S.; Qu, J.H. Hyperspectral pansharpening based on intrinsic image decomposition and weighted least squares filter. Remote Sens.
**2018**, 10, 445. [Google Scholar] [CrossRef] - Mookambiga, A.; Gomathi, V. Comprehensive review on fusion techniques for spatial information enhancement in hyperspectral imagery. Multidimens. Syst. Signal Process.
**2016**, 27, 863–889. [Google Scholar] [CrossRef] - Alparone, L.; Wald, L.; Chanussot, J.; Thomas, C.; Gamba, P.; Bruce, L. Comparison of pansharpening algorithms: Outcome of the 2006 GRS-S data-fusion contest. IEEE Trans. Geosci. Remote Sens.
**2007**, 45, 3012–3021. [Google Scholar] [CrossRef] - Vivone, G.; Alpharone, L.; Chanussot, J.; Mura, M.D.; Garzelli, A.; Licciardi, G.A.; Restaino, R.; Wald, L. A critical comparison among pansharpening algorithms. IEEE Trans. Geosci. Remote Sens.
**2015**, 53, 2565–2586. [Google Scholar] [CrossRef] - Du, Q.; Younan, N.H.; King, R.; Shah, V.P. On the performance evaluation of pan-sharpening techniques. IEEE Geosci. Remote Sens. Lett.
**2007**, 4, 518–522. [Google Scholar] [CrossRef] - Wald, L.; Ranchin, T.; Mangolini, M. Fusion of satellite images of different spatial resolutions: Assessing the quality of resulting images. Photogramm. Eng. Remote Sens.
**1997**, 63, 691–699. [Google Scholar] - Selva, M.; Santurri, L.; Baronti, S. On the Use of the Expanded Image in Quality Assessment of Pansharpened Images. IEEE Geosci. Remote Sens. Lett.
**2018**, 15, 1–5. [Google Scholar] [CrossRef]

Dataset | Reference HS Size | Simulated Size | Reference HS Spatial Resolution (SR) | Simulated SR | Band Number | Spectral Range |
---|---|---|---|---|---|---|

Washington DC | 200 × 200 | PAN 200 × 200 | 3m | PAN 3m | 191 | 0.4–2.5 $\mathsf{\mu}\mathrm{m}$ |

HS 50 × 50 | HS 12m | |||||

Salinas scene | 200 × 200 | PAN 200 × 200 | 3.7 m | PAN 3.7 m | 204 | 0.4–2.4 $\mathsf{\mu}\mathrm{m}$ |

HS 50 × 50 | HS 18.5 m | |||||

Dataset | Real HS Size | Real PAN size | Real HS SR | Real PAN SR | Band Number | Spectral Range |

Hyperion | 100 × 100 | 300 × 300 | 30 m | 10 m | 174 | 0.4–2.5 $\mathsf{\mu}\mathrm{m}$ |

Index | No Denoising | Average | Gaussian | Open | Closed | Open-Closing |
---|---|---|---|---|---|---|

RMSE | 0.0125 | 0.0119 | 0.0123 | 0.0121 | 0.0123 | 0.0112 |

SAM | 6.8506 | 6.8507 | 6.8506 | 6.8507 | 6.8506 | 6.8506 |

CC | 0.9156 | 0.9176 | 0.9176 | 0.9175 | 0.9176 | 0.9176 |

ERGAS | 25.5071 | 25.4812 | 25.4822 | 25.5052 | 25.4880 | 25.4804 |

Dataset | Index | Method | |||||
---|---|---|---|---|---|---|---|

GS | GFPCA | CNMF | HySure | Bayesian | HFWT | ||

Washington DC | RMSE | 0.0114 | 0.0139 | 0.0116 | 0.0097 | 0.0117 | 0.0112 |

SAM | 7.1069 | 9.8467 | 7.6438 | 6.8601 | 7.2586 | 6.8506 | |

CC | 0.8856 | 0.8179 | 0.8879 | 0.9018 | 0.8954 | 0.9176 | |

ERGAS | 36.0732 | 37.7303 | 34.3126 | 26.5678 | 29.3497 | 25.4804 | |

Salinas scene | RMSE | 0.0426 | 0.2224 | 0.0162 | 0.0163 | 0.0167 | 0.0138 |

SAM | 3.7807 | 2.9814 | 1.7586 | 1.7039 | 1.8015 | 1.5460 | |

CC | 0.8542 | 0.9429 | 0.9544 | 0.9583 | 0.9515 | 0.9625 | |

ERGAS | 4.3301 | 3.0843 | 2.6346 | 2.4451 | 2.9079 | 2.5312 |

Index | Method | |||||
---|---|---|---|---|---|---|

GS | GFPCA | CNMF | HySure | Bayesian | HFWT | |

RMSE | 0.0426 | 0.0451 | 0.0368 | 0.0398 | 0.0387 | 0.0358 |

SAM | 11.1037 | 15.6849 | 12.2755 | 12.8308 | 12.9605 | 9.1809 |

CC | 0.9296 | 0.9241 | 0.9631 | 0.9443 | 0.9504 | 0.9821 |

ERGAS | 15.3536 | 16.3011 | 11.1116 | 12.3441 | 12.2353 | 11.1109 |

Dataset | Method | |||||
---|---|---|---|---|---|---|

GS | GFPCA | CNMF | HySure | Bayesian | HFWT | |

Washington DC | 1.1764 | 2.3455 | 8.8369 | 43.3926 | 70.5347 | 6.9068 |

Salinas scene | 2.3953 | 4.8471 | 8.9498 | 55.4224 | 71.4972 | 7.1413 |

Hyperion | 2.6618 | 7.5118 | 23.3787 | 117.1729 | 158.7151 | 10.6448 |

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).