A COMPREHENSIVE REVIEW OF PANSHARPENING ALGORITHMS FOR GÖKTÜRK-2 SATELLITE IMAGES

In this paper, a comprehensive review and performance evaluation of pansharpening algorithms for GÖKTÜRK-2 images is presented. GÖKTÜRK-2 is the first high resolution remote sensing satellite of Turkey which was designed and built in Turkey, by The Ministry of Defence, TUBITAK-UZAY and Turkish Aerospace Industry (TUSAŞ) collectively. GÖKTÜRK-2 was launched at 18th. December 2012 in Jinguan, China and provides 2.5 meter panchromatic (PAN) and 5 meter multispectral (MS) spatial resolution satellite images. In this study, a large number of pansharpening algorithms are implemented and evaluated for performance on multiple GÖKTÜRK-2 satellite images. Quality assessments are conducted both qualitatively through visual results and quantitatively using Root Mean Square Error (RMSE), Correlation Coefficient (CC), Spectral Angle Mapper (SAM), Erreur Relative Globale Adimensionnelle de Synthése (ERGAS), Peak Signal to Noise Ratio (PSNR), Structural Similarity Index (SSIM) and Universal Image Quality Index (UIQI).


INTRODUCTION
Image fusion is an important research and application field for remote sensing imagery.In a general framework, image fusion can be considered as the extraction of high frequency details from a high spatial resolution image, and the injection of these components into a low spatial resolution image in order to obtain a high spatial and high spectral resolution image.If the high spatial resolution image is a PAN image, then this process is called 'pansharpening'.
A number of pansharpening review papers can be found in the literature.State-of-the-art and advanced pansharpening methods are critically reviewed for QuickBird satellite images, and context based methods are stated to provide the best result (Garzelli, 2004).Thomas et al. have showed that the Amélioration de la Résolution Spatiale par Injection de Structures (ARSIS) concept prevents spectral distortion of pansharpened images (Thomas, 2008).Pixel level image fusion is grouped into three categories, i.e., component substitution (CS) techniques, modulation based (MB) techniques and multiresolution analysis (MRA) based techniques in terms of the fusion mechanism they use by Jinghui et al. (Jinghui, 2010).Amro et al. have classified pansharpening techniques in terms of their utilization mechanism (Amro, 2011).Pansharpening methods included in commercially available software packages such as ERDAS, ENVI and ESRI are reviewed and evaluated (Zhang, 2012).Seventeen state-of-the-art and advanced methods for multispectral pansharpening are evaluated on two satellite images from IKONOS-2 and WorldView-2 and modulation transfer function (MTF) generalized Laplacian Pyramids (GLP) with High Pass Modulation (HPM) injection model (MTF-GLP-HPM) and Band Dependent Spatial Detail (BDSD) approach provide the best performances (Vivone, 2014).Comparisons of eleven different pansharpening approaches are conducted, and Bayesian sparse method is stated to provide the best performance (Loncan, 2015).Twenty one pansharpening techniques are presented and evaluated for the VNIR and SWIR bands of Sentinel-2, and it is obtained that Generalized Laplacian Pyramid * Corresponding author with MTF-matched filter and Context-Based Decision injection scheme (MTF-GLP-CBD) provides the best pansharpening result (Vaiopoulos, 2016).
This study is focused on pansharpening for GÖKTÜRK-2 satellite images.After BİLSAT and RASAT, our country's third remote sensing satellite GÖKTÜRK-2 is observing and monitoring the earth surface.GÖKTÜRK-2 was designed and produced in Turkey and provides a PAN image with 2,5m.spatial resolution, and four spectral bands (B, G, R, NIR) in 5m.spatial resolution each.There are various works in the literature for GÖKTÜRK-2 images.The results of orthorectification process of images from RASAT and GÖKTÜRK-2 satellites are presented with the aid of ERDAS LPS software in (Küpçü. 2014).The orbit tests of GÖKTÜRK-2 satellite system are provided in (Atak, 2015).Satellite image pre-processing steps such as ortorectification, radiometric calibration, geometric correction, pansharpening, contrast correction, MTF sharpening, and band registration for GÖKTÜRK-2 images are presented in (Teke. 2015) and these steps are demonstrated as a workflow in (Teke, 2016a).Detailed radiometric calibration research analysis is implemented for GÖKTÜRK-2 satellite images in (Teke, 2016b).GÖKTÜRK-2 images from the Zonguldak area are investigated for quality assessment in terms of geometric and radiometric features in (Topan, 2016).The process of absolute radiometric calibration of the GÖKTÜRK-2 satellite sensor was performed by using ground-based measurements for Tuz Gölü region in (Sakarya, 2016).The use of GÖKTÜRK-2 images for city planning, forestry, agriculture, land cover classification are presented in (Gürcan, 2016) and vegetation discrimination application is presented in (Kalkan, 2015).
There are also some pansharpening reviews in literature on GÖKTÜRK-2 satellite images.Intensity-Hue-Saturation (IHS), Brovey transform (BT) and Principal Component Analysis (PCA) pansharpening algorithms are analysed for RASAT and GÖKTÜRK-2 satellite images, and it is stated that for the former, PCA, and for the latter, BT gives the best performance (Ozendi, 2015).
This study also presents a review on pansharpening methods' performance for GÖKTÜRK-2 satellite images.However, the main distinction of this work with respect to other reviews on pansharpening for GÖKTÜRK-2 images is the number of pansharpening methods and the evaluation criteria utilized.This work is much more comprehensive in regards to previous work in both regards, and multiple GÖKTÜRK-2 images are evaluated both qualitatively and quantitatively to better highlight the performance variations.
The remainder of this paper is organized as follows.In Section two, the pansharpening approaches utilized in this work are presented.Experimental results are provided in Section three.The paper is concluded with some potential future studies in Section four.

PANSHARPENING APPROACHES
Brief descriptions of the pansharpening techniques utilized in this study are provided below.

Component Substitution (CS) based Pansharpening Algorithms
CS based pansharpening methods are based on the projection of the MS image into another space, replacing the component with the most spatial information with the histogram matched PAN image, and using inverse transformation to obtain the pansharpened image in the original space.
Because of the two stage (forward-backward) transformation matrix calculation, the classical CS fusion method generates a relatively high computational cost.The "general" CS fusion technique overcomes this shortcoming by using a simple linear equation proposed in (Dou, 2007): Although CS based approaches are fast for computing, easy for implementing and provide good spatial results, their main drawback is the spectral distortions they incur.
Intensity Hue Saturation (IHS) method is a standard CS method which generates the intensity component by taking the average of all MS bands, and then subtracts the PAN image, which is histogram matched to this I component, to obtain the spatial detail matrix.These details are injected into each MS band separately in order to obtain the pansharpened image (Te-Ming, 2001).
In BT technique, a synthetic PAN image is generated by taking the average of MS bands.By dividing the PAN image with this synthetic PAN image and multiplying the result by each MS band, the pansharpened image is acquired (Alan, 1987).
PCA technique makes the transformation of the MS image to feature space in order to obtain principal components.First principal component is assumed to contain most of the energy or most of the spatial information.This term is replaced by the PAN image which is histogram matched to this component.By using inverse PCA transform the pansharpened image is obtained (Pat, 1991).PCA transformation can be expressed in the general CS format (Te-Ming, 2001).
In Gram Schmidt (GS) based pansharpening technique, synthetic PAN image is acquired by using a GS mode, of which there are three available.In the first mode, the average of MS bands is taken as a synthetic image.In mode 2, low-pass filtered version of PAN image is received as a synthetic PAN image.For the last mode, a synthetic PAN image is achieved by using least square regression analysis (GS-LS).By subtracting this synthetic PAN image from the original PAN image, spatial details are obtained.Extracted details are injected into the MS bands which are upsampled to PAN resolution, in order to get the pansharpened image.The injection gain factor k g in GS pansharpening is defined as follows (Laben, 2000): where cov(.) is the covariance matrix and var(.) is the variance value.
In Hyperspherical Color Space (HCS) approach, MS image is first transformed into n-dimensional hyperspherical color space, and then, histogram matching is performed between the square of PAN image and the square of the intensity image.The histogram matched PAN image is replaced with the intensity component and the pansharpened image is obtained after inverse HCS transformation (Chris, 2010).HCS pansharpening approach can be easily described in terms of CS based general pansharpening model (Tu., 2012).
HPF pansharpening approach extracts the high frequency details from the PAN image using a HPF.As the next step, the details are multiplied with a weight factor, and these results are injected into MS bands separately for obtaining the pansharpened image (Ute, 2008).
University of New Brunswick (UNB) pansharpening technique uses a regression analysis between PAN and MS bands.In this way, a synthetic PAN image is obtained using the weighted summation of the MS bands.For obtaining the pansharpened image, PAN image is divided with the synthetic PAN image, and then multiplied by each MS band (Zhang, 1999).
Partial Replacement Adaptive Component Substitution (PRACS) technique firstly creates high and low resolution synthetic component images by using partial replacement in terms of linear regression and correlation coefficients among PAN and MS images.Thereafter, some coefficients are evaluated in order to minimize the distinction between PAN and MS image bands new adaptive CS-based fusion technique is described as follows (Choi, 2011): where k w is the statistical ratio DN (digital numbers) values, k  is the difference between the high and low resolution synthetic component image and the mean value difference between the mean of high and low resolution synthetic component image and k l L _ is the local instability regularization parameter.

Multiresolution Analysis (MRA) based Pansharpening Algorithms
CS based pansharpening approaches are good in spatial resolution, however, spectral distortion is their major drawback.MRA pansharpening techniques have been generated in order to overcome this drawback (Vivone, 2014).Some MRA based pansharpening approaches are introduced in the following part of this section.MRA based pansharpening can be formulated as follows (Vivone, 2015): where Smoothing Filter Based Intensity Modulation (SFIM) is a pansharpening technique that utilizes the ratio between PAN and its low pass filtered version.This ratio is then multiplied by each MS band separately in order to obtain the pansharpened image (Liu, 2000).
Wavelet Transform (WT) is a transformation used for decomposing an image into its high and low frequency components.In WT based pansharpening, first histogram matching is performed between PAN and each MS band, and then, WT is applied to the histogram matched PAN image.From the obtained decomposition form, only the approximation part is taken, and the other parts are set to zero.By using inverse WT, low resolution PAN image is obtained.Subtraction is applied from the original PAN image to create the detail matrix.Pansharpened bands are finally obtained by adding this matrix to each MS band (Fionn, 1995).
Additive Wavelet Luminance (AWL) is an MRA pansharpening approach originally proposed for three bands (RGB) MS images (Nuñez, 1999).AWL approach injects the high resolution details to the luminance component of the MS image.For this purpose, AWL technique first transforms the RGB image into IHS color space.After histogram matching between the PAN image and the intensity component, PAN and intensity component are decomposed into wavelet planes.And then, high spatial details from PAN image is added to intensity component's wavelet coefficients in order to obtain a merged intensity component.Finally, by using inverse IHS transform, pansharpened image is achieved.The original AWL approach can be easily extended to the more than three bands by using below equation (Otazu, 2005); where k MS ~ is the low resolution MS band upsampled to PAN image, and each MS band is divided to the mean of all MS bands.The multiplication of this proportional value tries to preserve spectral signature property.This technique is called proportional AWL (AWLP).
A Tróus Wavelet Transform (ATWT) is a multiresolution stationary wavelet decomposition algorithm.A tróus WT pansharpening is similar in operation to pansharpening by regular WT, with the exception of the utilized WT.A tróus WT is an undecimated wavelet transform and its low pass filter kernel is defined as follows (Nuñez, 1999): Burt and Adelson first developed the Laplacian Pyramid (LP) which is derived from the Gaussian Pyramid.The levels of Gaussian Pyramid constitute low-pass filtered version of an image, respectively.LP consist of detail images which are extracted from the subtraction of the same level image and its interpolated low pass version.Enhanced LP (ELP) have been introduced by the adding features that enforce zeroth level entropy and the lower correlation (Stefano, 1994).Generalized LP (GLP) is derived from ELP by adding the property of fractional scale number as a ratio (Bruno, 2012).In pansharpening by GLP, the pansharpened image is obtained after LP detail extraction from the PAN image, and then injecting details to the MS image upsampled to PAN resolution.

Variational Methods
P+XS technique uses the assumption that the geometry of spectral channels of MS image is related to the topographic map of corresponding PAN image.P+XS is a variational method and tries to minimize energy function which comprised of three components by using Gradient Descent algorithm.This variational framework carries out by minimizing the sum of integrals MS image and its low-pass filtered version and tangent vector multiplied gradient of MS bands and also the integral of the sum of subtraction the PAN image from alpha values multiplied MS bands.The reader is referred for detailed explanation about P+XS to (Baester, 2006).

Experimental Data Sets
The first image utilized in this work was acquired over İstanbul, Turkey, on 26.02.2015.A subset of 512 × 512 pixels is utilized for the experimental results.The PAN and RGB (the first three bands of MS) images of this dataset are provided in Fig. 1.It can be observed that this image has mixed urban/rural properties with lots of high frequency details.The second dataset was acquired over Büyükada, İstanbul on again the same date 26.02.2015.A subset of 300 × 300 pixels is utilized for the experimental results.The PAN and RGB (the first three bands of MS) images of this dataset are provided in Fig. 2. The scene includes sealand transition and urban features.The third data set was acquired over Alanya Tepebaşı on 12.11.2014.A subset of 300 × 300 pixels is utilized for the experimental results.The PAN and RGB (the first three bands of MS) images of this dataset are provided in Fig. 3.This scene includes an urban area and mountainous regions.For the experimental results, the original 4 band MS and PAN images acquired by GÖKTÜRK-2 have to be degraded in spatial resolution.For this purpose, both the MS and the PAN images are passed through Gaussian low-pass filters, and then decimated, i.e. downsampled, in spatial domain according to the resolution ratio between the original PAN image and MS bands, namely 2. Note that for the pansharpening methods to operate, the MS image to be sharpened has to be upsampled to the PAN resolution.For this purpose, the synthetsized MS image is upsampled to the synthetsized PAN image's size by using interp23tap filter.
The pansharpening methods presented in Section 2 are applied to the synthesized MS and PAN images, and the pansharpened MS images are compared with the original MS bands of the GÖKTÜRK-2 image by using a variety of performance metrics.
For RMSE, SAM and ERGAS, a smaller value indicates a better performance, whereas for the CC, PSNR, SSIM and UIQI metrics, the larger value indicates the better performance.
It should be noted that a better way to synthesize the images would be to directly follow the Wald's protocol such that the lowpass filters are characterized by the MTFs and spectral response functions (SRFs) of the sensor used for the acquisition, i.e.GÖKTÜRK-2 in this case.However, because the authors were not able to acquire the MTF and SRF values from the data providers, this protocol could not be used in this work.The lack of MTF and SRF values also prevent the use of some more recent pansharpening approaches, such as MTF-GLP-HPM, which depend on those values for operation.

Experimental Results
Quantitative experimental results are presented in Tables 1-3, whereas the visual results are provided in Figures 4-6, for the three utilized datasets, in order.
It can be observed from Table I, that overall, MRA pansharpening approaches have provided better performances for this dataset.The best performance has been obtained by the SFIM method, followed closely by ATWT based pansharpening, which in turn provides the best performance in terms of SAM.The worst pansharpening performances have been obtained by the CS methods of HCS, PCA and BT for this dataset.The PRACS method in the CS category provide the relatively best performance in its category.
Table II points out that for a more challenging dataset, the performances of all methods degrade.Whereas the worst performance are obtained when HCS and PCA again, followed closely by GS (in mode 1 only), the most powerful method in the CS category, namely PRACS has very slightly outperformed the MRA based SFIM method for most metrics.However, SFIM still outperforms PRACS for the ERGAS metric, whereas both methods are outperformed by ATWT when the PSNR metric is utilized.
For the Alanya dataset, it is observed from Table III that the situation for this dataset is very similar to the Büyükada dataset, in that PRACS is the best performing method, except for SFIM with ERGAS and ATWT with PSNR.The basic CS method are again outperformed significantly.For all datasets, the variational P+XS method has a medium performance in that it outperforms lots of methods while being outperformed in a similar manner.Visual results are provided for the Istanbul and Alanya datasets in Figure 4 and Figure 5, respectively.The third set of visual results could not be provided due to the page limit.From the qualitative results in Figures 4-5, it can be observed that while the visual outputs are somewhat similar for all methods, PRACS, SFIM, ATWT and AWLP methods provide clearer and sharped results, whereas the outputs of methods such as IHS, BT, PCA, GS1 and HCS are often blurred or have very slight color distortions.

CONCLUSIONS
In this letter, we present a comprehensive analysis of pansharpening algorithms for GÖKTÜRK-2 satellite images.It is performed visually and quantitatively by using statistical parameters and quantitative indexes.
In general, it can be stated that the CS based PRACS method and the MRA based SFIM method provides the best pansharpening performances for the utilized GÖKTÜRK-2 satellite images, followed closely by the MRA based ATWT method.The worst overall performances are obtained by the PCA, BT and HCS pansharpening approaches, which due to their simplicity and easy implementation, are among the most commonly utilized pansharpening approaches for RASAT and GÖKTÜRK-2 images in the literature.
These results are of course dependent on the acquisition sensor, and the results would potentially vary for other acquisition sensors, as the number of bands, the spectral ranges of the MS bands, and the spectral overlap between the PAN image and the MS bands have a direct effect on the pansharpening.
Another interesting point to note is that some simple CS methods have outperformed some of the more recent MRA based methods for two of the utilized datasets.The authors believe this may be caused because the spatial is most probably more effective for the four bands GÖKTÜRK-2 images over the spectral preservation, unlike some other sensors such as Landsat with a larger number of bands.Another possibility is that the MTFs have not been utilized in the synthetization process, which may result in a loss of an advantage for the MRA based methodologies.
Future studies may include taking GÖKTÜRK-2 MTF and SRF values, which the authors were not able to acquire for this work, into account for the experimental process, or utilizing pansharpening directly on the original acquired images and using without-reference quality metrics for performance assessment.
-sharpened image, P is the PAN image, k MS ~ is the MS image upsampled to the resolution of PAN, k g is the gain factor, k w is the modulation coefficient for the nth band, I is the intensity component of the MS image derived using the transformation and  is the spatial detail matrix derived using the PAN image.
image upsampled to the resolution of PAN image,  denotes element-wise multiplication, P is the PAN image and L P denotes that low pass version of the PAN image P .