Pansharpening Using Guided Filtering to Improve the Spatial Clarity of VHR Satellite Imagery

Abstract

Pansharpening algorithms are designed to enhance the spatial resolution of multispectral images using panchromatic images with high spatial resolutions. Panchromatic and multispectral images acquired from very high resolution (VHR) satellite sensors used as input data in the pansharpening process are characterized by spatial dissimilarities due to differences in their spectral/spatial characteristics and time lags between panchromatic and multispectral sensors. In this manuscript, a new pansharpening framework is proposed to improve the spatial clarity of VHR satellite imagery. This algorithm aims to remove the spatial dissimilarity between panchromatic and multispectral images using guided filtering (GF) and to generate the optimal local injection gains for pansharpening. First, we generate optimal multispectral images with spatial characteristics similar to those of panchromatic images using GF. Then, multiresolution analysis (MRA)-based pansharpening is applied using normalized difference vegetation index (NDVI)-based optimal injection gains and spatial details obtained through GF. The algorithm is applied to Korea multipurpose satellite (KOMPSAT)-3/3A satellite sensor data, and the experimental results show that the pansharpened images obtained with the proposed algorithm exhibit a superior spatial quality and preserve spectral information better than those based on existing algorithms.

1. Introduction

Very high resolution (VHR) satellite sensors, such as WorldView-3, Pléiades, and the Korea multipurpose satellite (KOMPSAT)-3/3A, provide panchromatic images with high spatial resolutions and multispectral images with low spatial resolutions. Generally, pansharpening is a methodology used to sharpen the spatial resolution or clarity of a multispectral image by adding spatial details from panchromatic images with high spatial resolutions [1]. Various pansharpening techniques have been proposed, following approximately two decades of research to extract spatial details from panchromatic images and then add those details through global/local methods [2,3,4,5]. An additional technique for enhancing image spatial resolution is hypersharpening, which is defined as enhancing the spatial resolution of a hyperspectral image by using multispectral or panchromatic image with high spatial resolution [6,7]. General pansharpening algorithms have been classified into component substitution (CS)-based and multiresolution analysis (MRA)-based methods depending on how the spatial details are generated [2]. CS-based algorithms generate pansharpened images by adding spatial details based on high-frequency information from panchromatic images with a high spatial resolution and synthetic intensity images with a low spatial resolution [8,9,10]. CS-based methods have the advantage of enhancing the spatial clarity of pansharpened images as the effects of aliasing, artifacts, and texture blurring are minimized during the pansharpening process [4]. The generalized intensity-hue-saturation (GIHS), Gram–Schmidt (GS), GS adaptive (GSA), and band-dependent spatial detail (BDSD) methods are representative CS-based pansharpening techniques [2,11,12]. Additionally, hybrid algorithms, such as partial replacement adaptive component substitution (PRACS) and generalized BDSD algorithms, have been developed in addition to various CS algorithms using global and local injection gains [13,14,15,16]. Alternatively, MRA-based pansharpening techniques generate pansharpened images using the differences in spatial characteristics between a panchromatic image with a high spatial resolution and a spatially degraded panchromatic image [2]. MRA-based algorithms, such as wavelet-based methods, high-pass filtering (HPF), generalized Laplacian pyramids with modulation transfer function (MTF)-matched filtering (MTF-GLP), and MTF-based algorithms using spatial principal component analysis (SPCA), efficiently preserve the spectral information of the original multispectral image [2,17,18,19,20,21]. However, some artifacts and texture blurring can occur in pansharpened images when applying MRA-based algorithms by utilizing the spatial dissimilarity between panchromatic and multispectral images [5,13]. Consequently, some studies have developed pansharpening algorithms to enhance the spatial clarity of pansharpened images based on MRA-based methods [22,23,24,25]. Nevertheless, most researchers have developed various pansharpening algorithms based on either CS or MRA aimed at generating multispectral images with a spatial resolution similar to that of a panchromatic image while preserving the spectral information of the former [2,18].

Various pansharpening algorithms have been proposed to solve the spectral distortion issues common among such techniques. Xu et al. [26] performed pansharpening to reduce the spectral distortion of pansharpened images by dividing panchromatic and multispectral images into several classes using the K-means algorithm and multiple regression equations. Restaino et al. [27] proposed a method for extracting synthetic panchromatic images by applying morphological operators in MRA-based fusion techniques and improving the spatial resolution compared to using traditional MRA-based techniques. Li et al. [28] proposed a segmentation-based pansharpening method for minimizing spectral distortion and increasing the sharpness of pansharpened images between vegetation and non-vegetation objects. Wang et al. [29] developed a new pansharpening model based on global and nonlocal spatial similarity regularizers to minimize local dissimilarities. Moreover, a pansharpening algorithm for preserving changes in vegetation cover was also proposed [30]. Furthermore, various injection gains using global, local, moving window, and segmentation methods have been applied to various pansharpening algorithms. Accordingly, a segmentation method was proposed by evaluating the time and accuracy associated with the calculation of injection gains [31]. Choi et al. [5] proposed a new hybrid pansharpening algorithm using local injection gains based on the normalized differential vegetation index (NDVI) to reduce computational costs.

Additionally, various pansharpening techniques based on deep learning techniques have been developed. Yang et al. [32] proposed PanNet, which is a deep learning architecture for solving the pansharpening problem associated with spectral and spatial preservation. Masi et al. [33] used a convolutional neural network composed of a three-layer architecture that includes several nonlinear spectral indices for pansharpening. Moreover, a learning method was developed for an efficient convolutional neural network by using a dilated multilevel block and deep residual network [34,35].

Recently, guided filtering (GF) was applied to generate spatial details and injection gains during the pansharpening process. In the improved adaptive intensity-hue-saturation (IAIHS) fusion algorithm, GF was used to compute the optimal weight of pansharpening [36]. Zheng et al. [37] utilized GF to properly add spatial details to imagery from the GaoFen-2 high-resolution imaging satellite. Liu and Liang [38] developed a pansharpening algorithm using GF to extract the missing spatial details of multispectral images by minimizing the difference between a panchromatic image and filtered output image. Additionally, GF based on three-layer decomposition was utilized for a pansharpening algorithm to efficiently extract spatial details from high-spatial resolution image [39]. In the abovementioned algorithms, multispectral images with a low spatial resolution are used as guidance images to optimize panchromatic images with a low spatial resolution; however, due to the differences associated with the time lag between panchromatic and multispectral image sensors, spatial dissimilarity could occur between panchromatic and multispectral images. Although these issues are important for improving the spatial clarity of pansharpened images, most methodologies have focused on minimizing spectral distortion rather than solving these problems.

Therefore, in this manuscript, we minimize the spatial dissimilarity between panchromatic and multispectral images and optimize the spatial clarity. In the proposed algorithm, GF is used to generate optimal multispectral images for pansharpening, in contrast to conventional GF-based pansharpening algorithms that use GF to extract spatial details. Additionally, the optimal panchromatic image possessing spatial characteristics similar to those of the multispectral image regenerated by GF is used to maximize the spatial details for pansharpening. Finally, during pansharpening, we modify the methodology for determining the NDVI-based optimal injection gains based on previous works [4,5]. In particular, the injection gains are optimized by a sigmoid function based on the characteristics of the NDVI, which exhibits a spectral pattern similar to general local injection gains, for KOMPSAT-3A. The proposed algorithm is then applied to satellite image products of KOMPSAT-3A to evaluate its performance on pansharpened products of full scenes. The new pansharpening algorithm based on GF and the modification of local injection gains are proposed in Section 2. In Section 3, the study area and materials are described. Section 4 and Section 5 provide an analysis and discussion of the experimental results based on a comparison of the quantitative and qualitative qualities of the pansharpened images obtained with our algorithm versus those obtained from existing state-of-the-art algorithms. Conclusions are presented in Section 6.

2. Guided Filtering (GF)-Based Pansharpening Algorithm

Recently, various studies of pansharpening algorithms using GF have been conducted. In this manuscript, we aim to generate an optimal multispectral image using GF, while most existing algorithms use GF to extract spatial details. The details of the proposed algorithm are shown in Figure 1.

2.1. Guided Filtering

Generally, GF algorithms are some of the most effective at removing noise from digital images and preserving edge information within images. A GF output image is obtained by a local linear model involving the filter output image $Q$ and guidance image $Y$ in a local window ${\omega }_k$, as shown in Equation (1) [40,41].

$$Q_i=a_k{Y}_i+b_k,i\in {\omega }_k$$

(1)

where $Q$ is modeled as the filter input image and the guidance image $Y$ by removing unwanted texture or noise. The linear coefficients of Equation (1) can be determined using the minimization of the squared difference $\mathrm{E}$ between the filter output image $Q$ and filter input image $X$ with a local ridge regularization parameter $\mathsf{\epsilon }$ in Equation (2).

$$\mathrm{E}\left(a_k,b_k\right)={{\displaystyle \sum }}_{i\in {\omega }_k}({(a_k{Y}_i+b_k-X_i)}^2+\epsilon a_k{}^2)$$

(2)

Therefore, the linear coefficients $a$ and $b$ of Equation (2) are determined by linear ridge regression according to Equations (3) and (4):

$$a_k=\frac{\frac{1}{\left|\omega \right|}{{\displaystyle \sum }}_{i\in w_k}{Y}_i{X}_i-{\mu }_k{\overline{X}}_k}{{\sigma }_k^2+\epsilon }$$

(3)

$$b_k={\overline{X}}_k-a_k{\mu }_k,{\overline{X}}_k=\frac{1}{\left|\omega \right|}{\displaystyle \sum }_{i\in {\omega }_k}{X}_i$$

(4)

where ${\mu }_k$ and ${\sigma }_k^2$ are the mean and variance of $Y$ in ${\omega }_k$, ${\overline{X}}_k$ is the mean of $X$ in ${\omega }_k$, and $\mathsf{\epsilon }$ is a GF regularization parameter [40]. After determining the GF coefficients using Equations (3) and (4), the filter output image can be defined based on Equation (5) through the reformulation of Equation (1):

$$Q_i=\frac{1}{\left|\omega \right|}{\displaystyle \sum }_{k:i\in {\omega }_k}\left(a_k{Y}_i+b_k\right)=\overline{a_i}{Y}_i+\overline{b_i}$$

(5)

where $\overline{a_i}=\frac{1}{\left|\omega \right|}{\displaystyle \sum }_{k\in {\omega }_i}{a}_k,\overline{b_i}=\frac{1}{\left|\omega \right|}{\displaystyle \sum }_{k\in {\omega }_i}{b}_k$, and $\left|\omega \right|$ is the number of pixels in the local window. For convenience, in this manuscript, the GF output image $Q$ is abbreviated following Equation (6) through the filter input image $X$, guidance image $Y$, window size $\omega$, and regularization parameter $\mathsf{\epsilon }$ [42]:

$$Q=G{F}_{\omega ,\epsilon }\left(X,Y\right)$$

(6)

2.2. GF-Based Pansharpening Algorithm

In general, the MRA-based pansharpened image ${\widehat{PS}}_k$ of the kth band can be determined using Equation (7):

$${\widehat{PS}}_k={\widetilde{MS}}_k+g_k\left(P-P^l\right),k=1,\dots ,N$$

(7)

where ${\widetilde{MS}}_k$ is the resized multispectral image of the kth band, $g_k$ denotes the injection gains of the pansharpening algorithm, and $P^l$ is a panchromatic image with a low spatial resolution. As noted in the previous chapter, CS-based pansharpening algorithms generate $P^l$ using linear combinations based on weight parameters and the relationship between the panchromatic and multispectral images. Alternatively, MRA-based algorithms use panchromatic images with a low spatial resolution by using various image degradation methods to obtain $P^l$.

2.2.1. Generation of an Optimal Multispectral Image with a Low Spatial Resolution Based on a Pansharpening Framework

In Equation (7), the spatial details are determined by the difference between $P$ and $P^l$, which has a low spatial resolution. In most GF-based pansharpening algorithms, the output filter image $Q_k$ of the kth band is applied as the synthetic intensity image of the pansharpening process to inject the optimal spatial details of the panchromatic image into each multispectral band [36,37,38]. These methods are similar to the MRA method insomuch as they generate the optimal panchromatic image with a low spatial resolution. However, some pansharpened images do not have abundant and clear spatial details, either because these details cannot be effectively injected into each band or because multispectral images with a lower spatial clarity than the original panchromatic image are employed. Figure 2a,b show the spatial characteristics of the target, which is composed of black and white tarps, acquired from panchromatic and resized KOMPSAT-3A multispectral images. In Figure 2b, a resized multispectral image is generated by cubic interpolation. As shown in Figure 2a,b, although the boundaries between the black and white tarps are very clear in the panchromatic image, various aliasing and noise sources around the boundaries of the tarps are included in the resized multispectral image. Although general pansharpening algorithms attempt to inject the spatial details of the panchromatic image into multispectral bands, these approaches also attempt to preserve the spectral information of the multispectral image as much as possible. Thus, the differences in the spatial details between the panchromatic and resized multispectral images can reduce the spatial clarity of pansharpened images. However, when GF is applied to a resized multispectral image, it is possible to generate an optimal multispectral image by removing these spatial characteristics. To generate filtered output images, previous studies of GF used multispectral images as guidance images and panchromatic images as input images. The GF technique uses a guidance image to generate an output image with noise-removed spectral characteristics similar to those of the input image. In this process, the output image has spatial characteristics similar to those of the guidance image. Therefore, when a multispectral image with a low spatial resolution is used as a guidance image, a panchromatic image with a low spatial resolution is generated. However, in this study, to remedy the degradation of the sharpness of the pansharpened image that occurs when the spatial characteristics of multispectral and panchromatic images are dissimilar, the original panchromatic image is used as a guidance image to generate multispectral images with spatial characteristics similar to those of the panchromatic image. Therefore, GF is applied to each band of the resized multispectral image using a panchromatic image as a guidance image. Then, the noise in the multispectral image regenerated by GF is removed using an MTF-matched filter. Figure 2c shows a multispectral image obtained by GF according to a resized multispectral image and a panchromatic image as a guidance image.

As shown in Figure 2c, the boundaries of the tarps and the edges of each object in the multispectral image generated by GF are clearer than those in the resized original multispectral image. This means that the optimal resized multispectral image for pansharpening can be generated using GF. When applying GF, the window size $\omega$ and regularization parameter $\mathsf{\epsilon }$ are set to 2 and 0.1, respectively, in reference to the results of Zheng et al. [37] related to pansharpening. In the proposed algorithm, pansharpening is performed by utilizing the filter output image obtained by GF as a multispectral image. Therefore, in the proposed algorithm, the general pansharpening framework in Equation (7) is revised to that in Equation (8):

$${\widehat{PS}}_k={\widetilde{MS}}_{GF,k}+g_k\left(P-P^l\right),k=1,\dots ,N$$

(8)

where ${\widetilde{MS}}_{GF,k}=G{F}_{2,0.1}\left({\widetilde{MS}}_k,P\right)$, which is the optimal resized multispectral image of the kth band obtained by GF. To generate ${\widetilde{MS}}_{GF,k}$, the original panchromatic image $P$ is used as the guidance image of GF, and the GF filter input image is set as the resized multispectral image ${\widetilde{MS}}_k$.

2.2.2. Local Injection Gains Based on a Sigmoid Function

The injection gain $g_k$ can also be an important factor in determining the pansharpening quality. Therefore, in this manuscript, the local injection gains are determined by modifying the NDVI-based local injection gains [5]. Xu et al. [26] indicated that the local injection gains of vegetated areas are different to those of non-vegetated areas. Additionally, Choi et al. [5] demonstrated that the spectral NDVI pattern reflects a high or moderate correlation to local injection gains based on moving windows. The local injection gain $g_k$ based on the NDVI is determined by Equations (9) and (10) [5]:

$$g_k={\left(-1\right)}^a\times NDVI+\overline{NDVI}+\sqrt{\frac{\sigma \left({\widetilde{\mathit{M}\mathit{S}}}_k\right)}{\sigma \left({\mathit{I}}_L\right)}\times {\left(C_k\right)}^3}$$

(9)

$$a=\{\begin{array}{c}1, \mathrm{if} \mathrm{corr}\left({\widetilde{MS}}_k, NDVI\right)<0\\ 0, \mathrm{if} \mathrm{corr}\left({\widetilde{MS}}_k, NDVI\right)>0\end{array}$$

(10)

where $\sigma \left(\mathrm{A}\right)$ is the standard deviation of $\mathrm{A}$ and $C_k$ is the high-frequency correlation value obtained by Laplacian filtering between ${\widetilde{MS}}_k$ and $I_L$. In Equation (9), $g_k$ is determined by the assumption that the NDVI and local injection gains are moderately or highly correlated. Therefore, the mean value of the NDVI is substituted as the global injection gain. However, in Equation (9), $C_k$, which is obtained by spatial correlation, might be underestimated due to spatial dissimilarity or misalignment between ${\widetilde{MS}}_{GF,k}$ and $I_L$. If $C_k$ is underestimated during pansharpening, the spatial details of the pansharpened image might be similar to those of the original multispectral image. Therefore, in this manuscript, we modify $C_k$, the maximum correlation value associated with spatial or spectral information, through Equations (11)–(13):

$$C_{max,k}=\max \left\{C_{spectral,k}, C_{spatial,k}\right\}$$

(11)

$$C_{spectral,k}=\frac{cov\left({\widetilde{MS}}_{GF,k},I_L \right)}{{\sigma }_{{\widetilde{MS}}_{GF,k}}{\sigma }_{I_L}}$$

(12)

$$C_{spatial,k}=\frac{cov\left(HPF({\widetilde{MS}}_{GF,k}),HPF\left(I_L\right) \right)}{{\sigma }_{HPF\left({\widetilde{MS}}_{GF,k}\right)}{\sigma }_{HPF\left(I_L\right)}}$$

(13)

where $cov\left(A,B\right)$ is the covariance between $A$ and $B,$ $HPF\left(A\right)$ is the high-pass-filtered image of $A$ using the Laplacian filter, and $I_L$ is the synthetic image generated by linear multiple regression between ${\widetilde{MS}}_{GF}$ and $P$.

Moreover, $g_k$ from Equation (9) might be underestimated if the NDVI values of the image have a large dynamic range. In such cases, some areas of $g_k$ have values that are negative or close to zero, and the spatial details may not be injected correctly based on Equation (8). Hence, overestimation could occur due to NDVI outliers or noise, causing the excessive injection of spatial details and the production of spectral distortions. Therefore, we reformulate Equation (6) using a sigmoid function, as described by Equations (14)–(16):

$$g_k=\left(\frac{1}{1+e^{-3\left\{{\left(-1\right)}^a\times NDVI+\overline{NDVI}\right\}}}+0.5\right)\ast g_G$$

(14)

$$g_G=\sqrt{\frac{\sigma \left({\widetilde{MS}}_{GF, k}\right)}{\sigma \left(I_L\right)}\times {\left(C_{max, k}\right)}^3}$$

(15)

$$a=\{\begin{array}{c}1, if corr\left({\widetilde{MS}}_{GF, k}, NDVI\right)<0\\ 0, if corr\left({\widetilde{MS}}_{GF, k}, NDVI\right)>0\end{array}$$

(16)

By using the sigmoid function in Equation (14), the spatial clarity in non-vegetated areas can be increased by setting a high $g_k$ value. In vegetated areas, we attempt to minimize spectral distortion by adjusting the value of $g_k$. Additionally, we minimize the underestimation of g in some regions by adjusting the parameters of the sigmoid function.

2.2.3. Extracting Spatial Details for Pansharpening

In this manuscript, we use a multispectral image based on GF. Therefore, the effect of GF is reflected when extracting spatial details from the original panchromatic image. Furthermore, by adjusting the panchromatic image, the problem of producing a pansharpened image that does not effectively increase the spatial clarity is avoided. First, $P^h$ with an increased spatial clarity is generated using $g_k$ and the MTF-matched filter. Equations (17) and (18) are used to generate $P^h$. In Equations (17) and (18), a constant value (0.5) is determined as the optimal value by experiments based on trial and error through 25 full scene products of KOMPSAT-3A:

$$P^h=P+\frac{1}{2}g{\prime }_k\left(P-P_{MTF}\right)$$

(17)

$$g{\prime }_k=\left(\frac{1}{1+e^{-3\left\{{\left(-1\right)}^{a }\times NDVI +\overline{ NDVI}\right\}}}+0.5\right)$$

(18)

where $P$ is the original panchromatic image and $P_{MTF}$ is the MTF-filtered image of $P$. In Equation (17), $P-P_{MTF}$ can be interpreted as the initial spatial details and $g_k$ is used to adjust the initial spatial details to ensure that excessive sharpening, which would result in spectral distortion during the pansharpening process, does not occur in vegetated areas. However, as $P_{MTF}$ from Equation (17) is only a filtered image and has not been subjected to the image upsampling process, we generated a synthetic panchromatic image, which can be used as $P^l$ in Equation (8), to extract the spatial details in the pansharpening process. Specifically, in this manuscript, since multispectral images were generated using GF, we performed pansharpening by generating synthetic panchromatic images $P_{GF}^l$ with characteristics similar to those of ${\widetilde{MS}}_k$. For this purpose, a low-resolution panchromatic image ${\tilde{P}}^l$ was generated through image downscaling and upsampling processes involving the MTF-filtered image $P_{MTF}$. Then, $P_{GF}^l$ was generated as a result of applying GF using ${\tilde{P}}^l$ as an input image and $P$ as a guidance image. As $P_{GF}^l$ is an image generated using the original panchromatic image, the proposed pansharpening algorithm can be classified as an MRA-based algorithm. However, through the extraction of spatial details by $P^h$ and $P_{GF}^l$, the spatial clarity of pansharpened images can be efficiently increased compared to that based on traditional CS- and MRA-based algorithms. In the case of non-vegetated areas, the spatial details obtained by the proposed method are clearer than those produced by traditional MRA-based algorithms, while the spatial details in vegetated areas produced by the existing and proposed algorithms are similar. Therefore, the spatial details extracted by the proposed algorithm can minimize spectral distortion in vegetated areas and effectively improve the spatial clarity. Finally, the proposed pansharpening algorithm in this paper can be defined as shown in Equation (19) by modifying Equation (8):

$$\begin{array}{cc}{\widehat{PS}}_k& ={\widetilde{MS}}_{GF,k}+g_k\left(P-P^l\right)\hfill \\ & ={\widetilde{MS}}_{GF,k}+g_k\left(P^h-P_{GF}^l\right)\hfill \\ & ={\widetilde{MS}}_{GF,k}+g_k\left(P-P_{GF}^l\right)+\frac{g_k^2}{2g_G}\left(P-P_{MTF}\right), k=1,\dots ,N\hfill \end{array}$$

(19)

3. Materials

The proposed pansharpening algorithm was applied to satellite imagery acquired by KOMPSAT-3A, which was launched by the Korea Aerospace Research Institute (KARI) on 26 May 2015. The specifications and characteristics of KOMPSAT-3A are described in

Table 1. The specifications of the Korea multipurpose satellite (KOMPSAT)-3A satellite sensor.

Sensor		KOMPSAT-3A
Multispectral resolution/size		2.2 m
Panchromatic resolution/size		0.55 m
Wavelength	Panchromatic	450–900 nm
	Blue	450–520 nm
	Green	520–600 nm
	Red	630–690 nm
	NIR	760–900 nm

.

In the experiment, two study areas, including targets based on tarps, were selected. The first site was located in the Salon region of France, and the second site was located in the Baotou region of China. The Salon region of France is complex, as it includes both urban and vegetated areas, whereas the Baotou region is composed primarily of cropland and natural terrain. In particular, all the satellite images used in the experiment were L1R products that were radiometrically corrected and covered full scenes.

Table 2. Descriptions of the experimental datasets.

	Site 1 (Salon)	Site 2 (Baotou)
Image size (panchromatic image)	$24,060\times 23,800$ (pixels)	$24,060\times 23,120$ (pixels)
Image size (multispectral image)	$6015\times 5950$ (pixels)	$6015\times 5780$ (pixels)
Acquisition date	16 July 2017	15 October 2017

describes the characteristics of the satellite images over the two sites, and Figure 3 illustrates the two study areas.

4. Experimental Results

4.1. Quality Assessment of Pansharpened Images

Various quality assessment methods and quality indices have been proposed to estimate the spectral and spatial quality of pansharpened images. In such cases, a multispectral image with a high spatial resolution that can be used as a reference should exist; however, such images are often unavailable. Therefore, many studies have proposed various quality assessment protocols to solve this problem. Such methods can be roughly divided into synthesis property, consistency property, and quality no reference (QNR) protocols [17,43]. To utilize the original multispectral image as a reference, the synthesis property method is used to generate a pansharpened image after downgrading the original multispectral and panchromatic images. Since the generated pansharpened image has the same spatial resolution as the original multispectral image, the pansharpened image and original multispectral image can be quantitatively compared based on the synthesis property. In the case of a consistency property protocol, the pansharpened image is generated using the original multispectral and panchromatic images, which are then spatially downgraded for comparison with the original multispectral image. The QNR protocol can be applied to a pansharpened image generated by the synthesis and consistency property methods, and the results are evaluated based on the relative similarity between the pansharpened and original multispectral images. Palsson et al. [17] showed that the consistency approach is the most reasonable evaluation method, although it has a similar tendency to the synthesis property. The QNR protocols could not efficiently reflect the spatial and spectral quality of pansharpened image. Therefore, this manuscript also performed evaluations using consistency and synthesis property.

4.2. Quality Indices for Estimating the Quality of a Pansharpened Image

To evaluate the spectral and spatial quality of pansharpened images using the consistency property approach, several spectral and spatial quality indices—namely, the Erreur Relative Globale Adimensionnelle de Synthèse (ERGAS), the spectral angle mapper (SAM), the universal image quality index (UIQI), and the correlation coefficient (CC) [2,17,44,45]—were employed to estimate the spectral quality. ERGAS estimates the global spectral/spatial error of pansharpened images using Equation (20) [2]:

$$\mathrm{ERGAS}=\frac{100}{\mathrm{R}}\sqrt{\frac{1}{\mathrm{N}}{\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{N}}{\left(\frac{\mathrm{RMSE}\left({\mathrm{I}}_{\mathrm{k}},{\mathrm{J}}_{\mathrm{k}}\right)}{\mathsf{\mu }\left({\mathrm{I}}_{\mathrm{k}}\right)}\right)}^2}$$

(20)

where ${\mathrm{I}}_{\mathrm{k}}$ is the reference multispectral image, ${\mathrm{J}}_{\mathrm{k}}$ is the pansharpened image, and $\mathrm{R}$ is the spatial resolution ratio between the multispectral and panchromatic images. The closer the ERGAS value is to zero, the less spectrally distorted the pansharpened image is. In the case of SAM, the average spectral difference in the angle between each pixel in the reference and pansharpened images is calculated using Equation (21) [2,17]:

$$\mathrm{SAM}=\arccos \left(\frac{\langle {\mathrm{I}}_{\left\{\mathrm{k}\right\}},{\mathrm{J}}_{\left\{\mathrm{k}\right\}}\rangle }{\Vert {\mathrm{I}}_{\left\{\mathrm{k}\right\}}\Vert \Vert {\mathrm{J}}_{\left\{\mathrm{k}\right\}}\Vert }\right)$$

(21)

where ${\mathrm{I}}_{\left\{\mathrm{k}\right\}}$ indicates a pixel vector of image $\mathrm{I}$ in the k-th band. Similar to ERGAS, the closer the SAM value is to zero, the less distorted the pansharpened image is. UIQI, which was developed by Wang and Bovik [44], reflects the loss of correlation, the luminance distortion and the contrast distortion and is calculated using Equation (22):

$$\mathrm{UIQI}=\frac{{\sigma }_{IJ}}{{\sigma }_I{\sigma }_J}\frac{2\overline{IJ}}{{\left(\overline{I}\right)}^2+{\left(J\right)}^2}\frac{2{\sigma }_I{\sigma }_J}{\left({\sigma }_I^2+{\sigma }_J^2\right)}$$

(22)

where ${\sigma }_{IJ}$ is the covariance of $I$ and $J$. The closer the UIQI value is to one, the less distorted the pansharpened image is. Additionally, CC is a representative spectral quality index for pansharpened images. Specifically, CC calculates the spectral similarity between the reference multispectral and pansharpened images using Pearson’s CC. The closer the CC value is to one, the greater the spectral similarity between the pansharpened image and the reference dataset [45].

4.3. Experimental Results and Analysis

To evaluate the performance of the proposed algorithm, we selected two pansharpening algorithms, namely, the GSA and MTF-GLP algorithms, for a comparison of the spectral and spatial quality of the pansharpened images [2,15]. In this manuscript, GFNDVI denotes the proposed GF-based pansharpening algorithm using local injection gains based on the NDVI. Figure 4 and Figure 5 show the pansharpening results according to each algorithm and detailed images of each pansharpened image. In the vegetated area (upper left area of Figure 4), the spectral distortion of the pansharpened images generated by the MTF-GLP is greater (Figure 4c) than that of the pansharpened images generated by the GSA and GFNDVI. Additionally, as shown in Figure 5, the colors of some cultivated areas in the GSA and MTF-GLP results are very bright; this brightness is caused by the excessive injection of spatial details into the vegetated area. Meanwhile, the images pansharpened by the GFNDVI have spectral and spatial characteristics similar to those of the original panchromatic and multispectral images. Furthermore, the pansharpened images by the GFNDVI have the best spatial clarity among the three techniques, as shown in Figure 4e and Figure 5e. This means that a multispectral image generated by GF can be utilized for pansharpening and that our methodology for extracting local injection gains is effective.

Table 3. Comparative pansharpening results corresponding to each region. ERGAS: erreur relative globale adimensionnelle de synthèse; SAM: the spectral angle mapper; CC: spatial correlation coefficient; UIQI: universal image quality index.

Region	Algorithm	Synthesis Property				Consistency Property
		ERGAS	SAM	CC	UIQI	ERGAS	SAM	CC	UIQI
Salon (France)	MTF-GLP	2.838	3.417	0.948	0.749	1.473	1.327	0.987	0.910
	GSA	2.630	3.352	0.955	0.777	1.453	1.315	0.987	0.945
	GFNDVI	2.728	3.850	0.951	0.742	1.290	1.644	0.991	0.936
Baotou (China)	MTF-GLP	0.826	0.884	0.966	0.867	0.434	0.493	0.993	0.952
	GSA	0.980	0.998	0.958	0.880	0.862	0.740	0.973	0.949
	GFNDVI	0.744	0.995	0.970	0.876	0.417	0.617	0.994	0.945

presents the results of the quality indices for the pansharpened images generated by each algorithm. For the calculation of the quality indices, in case of the consistency property, the downgraded pansharpening images were compared with the original multispectral images. In the case of the synthesis property, the pansharpened images were generated using downgraded panchromatic and multispectral images, and the evaluation index was calculated by comparing them with the original multispectral image. Therefore, the pansharpened image generated to apply the synthesis property should have the same spatial and spectral characteristics as the original multispectral and degraded panchromatic image. The pansharpened image used for the consistency property should be the same spatial characteristics as the original panchromatic image, while the degraded pansharpened image should be similar to original multispectral image. As shown in Table 3, the proposed method shows the best CC and ERGAS values, except for the results by synthesis property in the Salon region. However, in the case of SAM and UIQI, the MTF-GLP and GSA show the best results, respectively. The reason for this discrepancy is that the consistency property is used to quantify the difference between the pixel value from the original multispectral image. However, the proposed method applies GF to the original multispectral image to generate a spatially optimized multispectral image converted into characteristics similar to those of the original panchromatic image. As the original multispectral image is used as the reference for both the evaluation method by synthesis and consistency properties, pansharpened images generated by the GFNDVI algorithm can appear as if the evaluation indices for quantitative estimation might be decreased. This is due to the fact that GFNDVI aims to generate a pansharpened image which has similar spatial characteristics to the original panchromatic image. Additionally, since spectral and spatial quality are generally in a tradeoff relationship, spectral distortion may occur in a pansharpened image having better spatial clarity. Therefore, the evaluation indices of the pansharpened images generated for the Salon region become low, since the spatial clarity is greatly improved, as shown in Figure 4e, and the spectral characteristics of spatially enhanced area are emphasized even more. Nevertheless, considering the best ERGAS and CC values for both the Salon and Baotou regions, the proposed method effectively preserves the spectral information of the original multispectral image.

Additionally, to quantitatively analyze the spatial clarity of the pansharpened images, we analyzed the spatial characteristics of edge targets existing within the image. The results of enlarging the edge target area existing in the pansharpened image are shown in Figure 6 and Figure 7; evidently, the pansharpened images generated by the existing techniques do not show an edge target with definite linearity. In particular, the MTF-GLP, which is an MRA-based pansharpening algorithm, greatly distorts the edge characteristics of the target. Aliasing and blurring are observed in the edges around the target, even in the case of the images fused by the GSA. However, the pansharpened image generated by the GFNDVI effectively represents the edge information of the target, as shown in Figure 6e and Figure 7e. In particular, as shown in the specific area of the red rectangle in Figure 6a and Figure 7a, the edge lines between the black and white areas in the pansharpened images generated by the MTF-GLP and GSA include the effects of aliasing and artifacts. However, the edge target in the pansharpened image produced by the proposed algorithm does not exhibit artifacts. Figure 8 presents magnified views of the area encompassed by the red rectangle in Figure 7a. As shown in Figure 8b,c, the images pansharpened by the MTF-GLP and GSA exhibit aliasing around the edges between the black and white areas of the edge target. Furthermore, the image pansharpened by the MTF-GLP displays spectral distortion around the cross line of the edge target, as shown in Figure 8b. However, the image pansharpened by the proposed algorithm has similar spatial characteristics to those of the original panchromatic image, as shown in Figure 8a,d. Therefore, our proposed algorithm preserves the spatial information of the original panchromatic image during the pansharpening process, while minimizing spectral distortion.

To quantitatively verify the spatial quality results of Figure 6 and Figure 7, the edge of each target was extracted; furthermore, the signal-to-noise ratio (SNR) and the Nyquist value of the MTF based on both an edge spread function (ESF) and a line spread function (LSF) were calculated [46,47]. The reference dataset for the SNR and Nyquist values was based on the original panchromatic image. Since the center of the target is composed of four edges, four edges are extracted from the image, after which the SNR and Nyquist values for the edges are calculated. If the edges in the image are clear, the SNR and MTF-Nyquist values calculated through those edges are high.

Table 4. Comparative spatial clarity results corresponding to each region. SNR: signal-to-noise ratio; MTF-Nyquist: Nyquist value based on the MTF.

Region	Algorithm	SNR (dB)	MTF-Nyquist (%)
Salon	Panchromatic	67.55	17.11
	MTF-GLP	52.17	14.74
	GSA	51.67	16.37
	GFNDVI	63.11	17.90
Baotou	Panchromatic	47.65	26.63
	MTF-GLP	42.21	17.79
	GSA	42.29	20.76
	GFNDVI	47.37	23.91

shows the average SNR and Nyquist values for each edge. As shown in Table 4, the SNR and Nyquist values for the edges in the pansharpened images generated by the GFNDVI are higher than those generated by the MTF-GLP and GSA. Table 4 also shows that the SNR and MTF-Nyquist values are similar to those of the original panchromatic image. Therefore, the existing pansharpening algorithms generate aliasing and blurring effects along edge boundaries; however, the proposed algorithm effectively reflects the spatial characteristics of the original panchromatic image.

5. Discussion

The experimental results confirmed that the proposed GFNDVI technique produces similar or superior pansharpened images in terms of the spectral and spatial quality in comparison with existing pansharpening techniques. However, the aim of this study is to generate spatially optimal pansharpened images by removing spatial dissimilarity from multispectral and panchromatic images. Therefore, in this section, the efficiency of the revised algorithm in extracting the local injection gains is discussed. The proposed approach for extracting the local injection gains derives the optimal variables that are neither overestimated nor underestimated. To evaluate these claims, we compared the extracted results of the local injection gains using Equation (9) with the results using the proposed method. The average, maximum, and minimum values of the local injection gains for each band extracted through each technique are shown in

Table 5. Statistical characteristics of local injection gains according to algorithms used to extract the local injection gains.

Algorithm	Band	Average	Max. Value	Min. Value
Method by Equation (9) [5]	Blue	0.6931	2.1160	0.1158
	Green	1.1217	2.5444	0.5444
	Red	1.2281	2.6077	0.6508
	NIR	1.4841	2.0614	0.0614
Proposed method by Equation (14)	Blue	0.6929	1.0301	0.4508
	Green	1.1213	1.6671	0.7295
	Red	1.2276	1.8252	0.7987
	NIR	1.4847	2.0030	0.7625

, which demonstrates that the averages of both methods are similar but that the local injection gains of the blue and NIR bands generated by Equation (9) are close to zero. Additionally, the maximum value is excessively large. However, we can confirm that this tendency is eliminated in the case of the proposed scheme using a sigmoid function.

Figure 9 presents a histogram plot of the local injection gains for the blue band for the Salon region of France shown in Figure 3a. As shown in Figure 9b, the histogram plot generated by the proposed method shows a decreased dynamic range due to the minimization of overestimated and underestimated values in comparison with Figure 9a. Therefore, the technique for extracting local injection gains proposed in this manuscript is slightly more stable than the existing technique.

Moreover, the spatial quality of images pansharpened using an ESF and LSF is further analyzed. Figure 10 and Figure 11 show the ESF and LSF along and across the panchromatic image in the blue band according to each algorithm in the Salon and Baotou regions, respectively. As shown in Figure 10 and Figure 11, the ESF curves of the pansharpened image generated by the GFNDVI present a pattern similar to that of the original panchromatic image compared with those of the pansharpened images generated by the GSA and MTF-GLP. Additionally, the LSF curves of the pansharpened images obtained by the GSA and MTF-GLP have a relatively wide full width at half maximum (FWHM), and the distortion of the LSF curve is large compared to that of the original panchromatic image. These results confirm that the linearity of the edge target generated through the proposed algorithm is associated with a small error and high SNR and MTF values. This trend was common at both experimental sites, as shown in Figure 10 and Figure 11. Therefore, the proposed GFNDVI algorithm effectively preserves the spatial clarity of original panchromatic images in the pansharpening process. Notably, the proposed method yields similar and ideal ESF and LSF curve shapes compared with the existing techniques. Furthermore, the proposed technique retrieves high SNR and MTF-Nyquist values.

6. Conclusions

In this manuscript, a new GF-based pansharpening algorithm is proposed to minimize spectral and spatial distortion in pansharpened images caused by spatial dissimilarities due to the differences between panchromatic and multispectral images related to the time lag between the sensors. Specifically, the proposed algorithm is focused on maintaining the spatial clarity of the original panchromatic image and minimizing spectral distortion within the pansharpened image. The main cause of a decrease in spatial clarity is that a resized multispectral image does not have the same spatial characteristics as a panchromatic image. Therefore, GF is used to generate an optimal multispectral image with the same spectral characteristics as the resized multispectral image and spatial characteristics similar to those of the panchromatic image. Additionally, to extract the local injection gains specific to the resized multispectral image based on GF, the existing injection gains were optimized using a sigmoid function. The quality of the pansharpened image generated through the proposed technique was analyzed based on the spectral and spatial characteristics of existing pansharpening image evaluation techniques and targets within images. The experimental results show that the proposed method yields less spectral distortion and better spatial clarity than conventional pansharpening algorithms. The computational costs of extracting local injection gains and the pansharpening model are similar to those of the general GSA and MTF-GLP algorithms; however, further works using parallel processing or graphics processing units will be needed since GF requires a relatively high computational cost.