An Approach for the Pan Sharpening of Very High Resolution Satellite Images Using a CIELab Color Based Component Substitution Algorithm

Recent very high spatial resolution (VHR) remote sensing satellites provide high spatial resolution panchromatic (Pan) images in addition to multispectral (MS) images. The pan sharpening process has a critical role in image processing tasks and geospatial information extraction from satellite images. In this research, CIELab color based component substitution Pan sharpening algorithm was proposed for Pan sharpening of the Pleiades VHR images. The proposed method was compared with the state-of-the-art Pan sharpening methods, such as IHS, EHLERS, NNDiffuse and GIHS. The selected study region included ten test sites, each of them representing complex landscapes with various land categories, to evaluate the performance of Pan sharpening methods in varying land surface characteristics. The spatial and spectral performance of the Pan sharpening methods were evaluated by eleven accuracy metrics and visual interpretation. The results of the evaluation indicated that proposed CIELab color-based method reached promising results and improved the spectral and spatial information preservation.


Introduction
The earth observation satellites with very high resolution (VHR) optical sensors provide a multispectral (MS) image and a panchromatic (Pan) image that are acquired simultaneously in order to provide essential accommodation between spectral and spatial resolution, which is an important consideration for optical satellite sensors due to their physical limitations [1,2]. Spectral diversity is important for modeling the spectral characteristics of different land cover/use classes and identifying them; on the other hand, spatial information is very crucial for identifying spatial details and geometric characteristics. The Pan image provides high spatial resolution with a single, wide range spectral band, whereas the MS image provides several spectral bands in different sections of the electromagnetic spectrum with low spatial resolution in order to meet the abovementioned requirements.
The fusion of Pan and MS images that are acquired over the same area from the single or multiple satellite system is referred to as Pan sharpening. The main aim of Pan sharpening is to create a high-resolution MS image, having the spatial resolution of Pan but preserving the spectral characteristics of MS [3]. Unlike the challenging problem of multi-sensor data fusion, single sensor Pan sharpening does not need image-to-image registration, as the Pan and MS sensors are mounted on the same platform and the images are acquired simultaneously with well-matching viewing geometry [4]. Several earth observation satellites, such as Geo-Eye, OrbView, QuickBird, WorldView, Pléiades and

State of Art
A brief review of the state of the art methods that used for comparative analysis with respect to the proposed method are presented in Table 1. Table 1. The brief review of state of the art Pan sharpening methods.

Method
Description References

IHS Method
In this system, the total amount of the brightness in one color is represented through intensity channel. The wavelength property of the color and the purity of the color are represented by hue and saturation respectively. [30,31] EHLERS (FFT-Enhanced IHS Transform) Method The fundamental idea of this method is modifying the Pan image in the way that it looks more similar to the intensity component of MS image. This method uses FFT (fast fourier transform) filtering for partial replacement instead of entire replacement of the intensity component. [32] NNDiffuse (Nearest-neighbor diffusion-based) Method This method, considers each pixel spectrum as a weighted linear combination of spectra of its sideward neighboring super pixels in the Pan sharpened image. Algorithm uses various factors like intensity smoothness (σ), spatial smoothness (σ_s) and pixel size ratio for conducting the Pan sharpening [33,34] GIHS (Generalized IHS) Method Directly applying of IHS method needs many multiplication and addition operations, which makes the Pan sharpening operation computationally inefficient. GIHS method develops a computationally efficient Pan sharpening method, which does not require coordinate transformation [35] Gram-Schmidt Method This method uses the Gram-Schmidt orthogonalization for converting the original low-resolution MS bands, which are linearly independent vectors, into a set of orthogonal vectors. The first vector in the orthogonal space is considered as simulated Pan image, which produced by weighted aggregation of the consecutive original MS bands. [21,22] Hyperspherical Color Sharpening Method The Hyperspherical Color Sharpening method (HCS) is a Pan sharpening method designed for WorldView-2 sensor imagery and can be applied to any MS data containing 3bands or more. HCS approach is based on transforming original color space to hyperspherical color space. [36,37]

Methodology and Accuracy Assessment
This research proposes a new Pan sharpening method that relies on the CIELab transform, which modifies the color components of images to be used for Pan sharpening of VHR satellite images. The flowchart of the proposed method is given in Figure 1, and the details of the proposed method are described in the sub sections.

Multispectral Image Transform to CIELab Color System
The RGB color system was designed in such a way that it includes nearly all primary colors and can be comprehended by human vision. Nevertheless, it is a tough task to deal with RGB color due to strong correlation between its components [29]. In this study, a uniform and complete color model, the CIELab color system, was used. In this uniform color system, a variation in the coordinates of the color component provides the same amount of variation in the luminance and saturation components [10]. Besides, this color space is projected to draw human vision, unlike the RGB and CMYK (cyan, magenta, yellow, black) color spaces [38]. The CIELab color system was used in the proposed Pan sharpening approach in order to reduce the spectral distortion, while maintaining the color perception of human vision [39]. Appl

Multispectral Image Transform to CIELab Color System
The RGB color system was designed in such a way that it includes nearly all primary colors and can be comprehended by human vision. Nevertheless, it is a tough task to deal with RGB color due to strong correlation between its components [29]. In this study, a uniform and complete color model, the CIELab color system, was used. In this uniform color system, a variation in the coordinates of the color component provides the same amount of variation in the luminance and saturation components [10]. Besides, this color space is projected to draw human vision, unlike the RGB and CMYK (cyan, magenta, yellow, black) color spaces [38]. The CIELab color system was used in the proposed Pan sharpening approach in order to reduce the spectral distortion, while maintaining the color perception of human vision [39].
The design of the CIELab color system is based on Hering's theory, which indicates that only red (R), green (G), blue (B) and yellow (Y) are unique among the thousands of colors that are used to characterize the hue component [40]. Although, other colors can be produced using these unique colors (for example it is possible to obtain orange by mixing red and yellow), they (R, G, B and Y) can be described only with their own name. R, G, B and Y, with black (B) and white (W), constitutes a color system with six basic color properties and three opponent pairs: R/G, Y/B and B/W. The opponency idea rises from observation upon colors attributes, which proves no color could be characterized using both blue and yellow or red and green together [41]. A blue shade of yellow does not exist. These three opponent pairs are represented in the form of a three-dimensional color space, as illustrated in Figure 2. In this figure, the vertical axis L* represents the luminance, in which perfect black is represented by 0 and perfect white is represented by 100. The a* and the b* are the axes that are perpendicular to luminance indicated chromaticity, and stand for redness/greenness and yellowness/blueness, respectively. Positive values represent redness (for a* component) and yellowness (for b* component), whereas greenness and blueness are denoted with negative values.  The design of the CIELab color system is based on Hering's theory, which indicates that only red (R), green (G), blue (B) and yellow (Y) are unique among the thousands of colors that are used to characterize the hue component [40]. Although, other colors can be produced using these unique colors (for example it is possible to obtain orange by mixing red and yellow), they (R, G, B and Y) can be described only with their own name. R, G, B and Y, with black (B) and white (W), constitutes a color system with six basic color properties and three opponent pairs: R/G, Y/B and B/W. The opponency idea rises from observation upon colors attributes, which proves no color could be characterized using both blue and yellow or red and green together [41]. A blue shade of yellow does not exist. These three opponent pairs are represented in the form of a three-dimensional color space, as illustrated in Figure 2. In this figure, the vertical axis L* represents the luminance, in which perfect black is represented by 0 and perfect white is represented by 100. The a* and the b* are the axes that are perpendicular to luminance indicated chromaticity, and stand for redness/greenness and yellowness/blueness, respectively. Positive values represent redness (for a* component) and yellowness (for b* component), whereas greenness and blueness are denoted with negative values.

Multispectral Image Transform to CIELab Color System
The RGB color system was designed in such a way that it includes nearly all primary colors and can be comprehended by human vision. Nevertheless, it is a tough task to deal with RGB color due to strong correlation between its components [29]. In this study, a uniform and complete color model, the CIELab color system, was used. In this uniform color system, a variation in the coordinates of the color component provides the same amount of variation in the luminance and saturation components [10]. Besides, this color space is projected to draw human vision, unlike the RGB and CMYK (cyan, magenta, yellow, black) color spaces [38]. The CIELab color system was used in the proposed Pan sharpening approach in order to reduce the spectral distortion, while maintaining the color perception of human vision [39].
The design of the CIELab color system is based on Hering's theory, which indicates that only red (R), green (G), blue (B) and yellow (Y) are unique among the thousands of colors that are used to characterize the hue component [40]. Although, other colors can be produced using these unique colors (for example it is possible to obtain orange by mixing red and yellow), they (R, G, B and Y) can be described only with their own name. R, G, B and Y, with black (B) and white (W), constitutes a color system with six basic color properties and three opponent pairs: R/G, Y/B and B/W. The opponency idea rises from observation upon colors attributes, which proves no color could be characterized using both blue and yellow or red and green together [41]. A blue shade of yellow does not exist. These three opponent pairs are represented in the form of a three-dimensional color space, as illustrated in Figure 2. In this figure, the vertical axis L* represents the luminance, in which perfect black is represented by 0 and perfect white is represented by 100. The a* and the b* are the axes that are perpendicular to luminance indicated chromaticity, and stand for redness/greenness and yellowness/blueness, respectively. Positive values represent redness (for a* component) and yellowness (for b* component), whereas greenness and blueness are denoted with negative values.  The L*, a* and b* values are computed using XYZ values. The XYZ system that was based on the RGB color space was presented by the International Commission on Illumination, CIE (Commission international de l'éclairage), in the 1920s and patented in 1931. The difference of RGB and XYZ lies in the light sources. The R, G and B elements are real light sources of known characteristics, whereas the X, Y and Z elements are three theoretical sources, which are selected in a way that all visible colors can be defined as a density of just-positive units of the three primary sources [10].
Occasionally, the colorimetric calculations with the use of color matching functions produce negative lobs. This problem can be solved by transforming the real light sources to these theoretical sources. In this color space, red, green and blue colors are more saturated than any spectral RGB. X, Y and Z components, represent red, green and blue colors respectively. RGB to XYZ and its reverse transformations can be performed by following equations: (2 Lab system is calculated by the following equations [23]: Where And where X n is the tristimulus value of a perfect white object color stimulus, which the light reflected from a perfect diffuser under the chosen illuminant. F Y and F Z values are calculated in the same way as F X .

Pan-Sharpening
In the Pan sharpening procedure, the MS image should be resampled to the same pixel size of the Pan image before converting it to CIELab color space. In this research, the bicubic interpolation method was used to resample 2 m resolution Pleiades MS images into 50 cm resolution to match the pixel size of Pleiades Pan image. This resampled dataset is used in all Pan sharpening methods used in this research, including the proposed one. After converting the MS image from RGB to CIELab space, the Pan sharpening process continues with replacing the Pan image with the L* component. Unlike the proposed method in [29] study, there is no need for color space conversion of Pan image in the proposed method, which leads to low computation and less data distortion. Before replacing the L* band of MS with the Pan image, there is a histogram matching step that could be considered as preprocessing step. After resampling the MS image to the same size of Pan image and converting the MS image color space, the histogram of Pan image has to be matched with the histogram of L* component in order to minimize the spectral differences [42]. For performing histogram matching task, mean and standard deviation normalizations were used [43]: where Pan HM stands for histogram matched Pan image, µ stands for mean and σ represents standad deviation. After these preprocessing steps, the L* component is replaced with a Pan image. The Pan sharpened image is then produced by implementing inverse conversion of CIELab color system on the Pan*a*b* image and results in a new MS image with high spatial resolution.

Accuracy Assessment
Several metrics were proposed to assess the accuracy of Pan-sharpened images that use the precise, high-resolution MS image as a reference image. In this research, the first seven metrics provided in Table 2 were used for the spectral quality assessment, while the later four metrics were used for spatial quality assessment of the results. Although metrics provide important quantitative insights about the algorithm performance, qualitative assessment of the color preservation quality and spatial improvements in object representation is required. Thus, results obtained from the Pan sharpening algorithms were also evaluated with visual inspection.
Lower (near to zero) [45] ERGAS Relative dimensionless global error synthesis (ERGAS) is used to calculate the accuracy of Pan sharpened image considering normalized average error of each band of the result image.
Lower (near to zero) [45] SAM Spectral Angle Mapper (SAM) represents the spectral similarity between the Pan sharpened and reference MS image using the average spectral angle Lower (near to zero) [36] RASE Relative average spectral error (RASE) is an error index to calculate average performance of Pan sharpening algorithm into spectral bands.
Lower (near to zero) [46] PSNR Peak signal-to-noise ratio is widely used metric for comparison of distorted (Pan sharpened) and original (reference) image. PSNR = 20 log 10 Higher Value [47] QAVG The Average Quality index based on quality index is used to model the difference between reference and Pan sharpened images as a combination of three different factors: loss of correlation, luminance distortions and contrast distortion. As QI can only be applied to one band, the average value of three or more bands (QAVG) is used for calculating a global spectral quality index for multiband images.
Higher Value (Close to 1) [48] SSIM Structural Similarity index (SSIM) is a method for measuring the structural similarity between reference and Pan sharpened images. This method compares the local patterns (luminance, contrast and structure) using means and standard deviations of two images.
Higher Value [2] CC To assess the spatial quality of Pan sharpened images, the correlation coefficient between the Pan image and the intensity component of the Pan sharpened image is used.
Higher Value (Close to 1) [49] ZHOU index Zhou's spatial index uses a high frequency Laplacian filter for extracting high frequency information from both Pan and Pan sharpened images. Correlation coefficient is then calculated between filtered Pan image and each band of Pan sharpened image. The average of calculated cc is considered as spatial quality index.
SRMSE Sobel based RMSE (SRMSE) is an index for spatial accuracy assessment that uses absolute edge magnitude difference of Pan and Pan sharpened image. This index utilizes 3 × 3 vertical and horizontal Sobel filter kernels for calculating the gradient of edge intensities. RMSE then calculated between Pan and Pan sharpened edge magnitude images.
Lower (near to zero) [51] Sp-ERGAS Spatial ERGAS (Sp-ERGAS) is an index for spatial quality assessment of Pan sharpened image, which uses spatial RMSE for assessment procedure Lower (near to zero) [52]

Dataset
The primary product of Pléiades satellite images were used for performing experimental analysis of the proposed method. The Pléiades program, which was launched by CNES (the French space agency), is the optical Earth imaging component of French-Italian ORFEO (Optical and Radar Federated Earth Observation). The Pléiades constellation consists of two satellites with VHR optical sensors. The Pléiades 1A launched on 17.12.2011 and Pléiades 1B launched on 2.12.2012. Both of the satellites provide 0.5 m spatial resolution for the Pan sensor and 2 m for MS sensor with 20 km swath width and 12 bit radiometric resolution [44].
The dataset used in this research consists of three Pleiades image scenes, which cover different landscape characteristics ( Table 3). The locations of scenes are provided in Figure 3. Ten different sub frames were selected from these image scenes, in order to evaluate the performance of Pan sharpening methods for varying landscape characteristics and seasonal conditions ( Table 4). Four sub frames cover rural areas that contain forests with different types of trees with varying heights and barren roads between them, and mountains and agricultural fields. Three sub frames cover urban areas that include complex buildings, roads and highways. Finally, three sub frames cover sub-urban areas that include all the complex buildings, factories, vegetation, roads, trees and bare soil parts to perform a precise survey on effects of the proposed method in Pan sharpening of vegetation, impervious and soil surfaces simultaneously.      Four sub frames cover rural areas that contain forests with different types of trees with varying heights and barren roads between them, and mountains and agricultural fields. Three sub frames cover urban areas that include complex buildings, roads and highways. Finally, three sub frames cover sub-urban areas that include all the complex buildings, factories, vegetation, roads, trees and bare soil parts to perform a precise survey on effects of the proposed method in Pan sharpening of vegetation, impervious and soil surfaces simultaneously.

Performing the Algorithm and Accuracy Assessment
To measure the performance of Pan sharpening results using the metrics that were presented in Section 3.3, the Wald protocol was used, due to lack of reference a high-resolution MS image [53]. According to Wald protocol, all Pan sharpening experiments were done using degraded datasets, which are produced by decreasing spatial resolution of the original dataset (reduce MS and Pan, respectively, to 8 m and 2 m). The Pan sharpening results obtained that way, can be compared with the original MS images for an accuracy assessment procedure. In this paper, six Pan sharpening methods and eleven accuracy indexes are evaluated to perform a comparative accuracy assessment of the proposed method. The numerical results of the accuracy indexes are presented in Tables A1-A6. The visuals belonging to Pan sharpening results of ten frames are presented in Figures A1-A10. In each figure, parts a and b are the original Pléiades MS and Pan images, respectively. Parts c, d and e are the Pan sharpened results from the CIElab, GIHS and GS methods, respectively. The Pan sharpened images from the HCS, IHS, NNDiffuse and Ehlers methods are shown in parts f-i, respectively.

Experimental Results from Rural Areas
Figures A1-A4 belongs the Pan sharpening results of the rural test sites (frame F1 from D1 dataset, frames F5 and F6 from D2 and Frame F10 from D3 datasets). Each figure belongs to a representative part from the whole image focusing on rural areas and presents visual comparison different Pan sharpening techniques.
The visual comparison of the Pan sharpening methods reveals that spatial resolution of MS images improved significantly in all methods. As for spectral information, parts c, e and h show that the CIELab GS and NNDiffuse methods protect the spectral characteristics better; specifically, for the bands belonging to the visible region. The color-based visual interpretation in vegetated and forest areas in Figures A2-A4 inform us that the Pan-sharpened and original MS images are very similar to each other for GS and CIELab methods. Similar comments can be made on NNDiffuse and CIELab methods in Figure A1. On the other hand, visual comparison of part a with parts d, f, g and i reveals that the remaining four methods were not able to preserve the spectral characteristics of vegetated and forest areas. Particularly, IHS, Ehlers and HCS methods inherited the high frequency impact over vegetated area and could not preserve original spectral/color information for the first test site. In addition, the result of the HCS method is more blurred than the others. The GIHS and-in some cases-the NNDiffuse methods, preserved the color information better than the IHS, Ehlers and HCS; nevertheless, observable spectral distortion is apparent in their resulting products. The GS method has good performance in the case of vegetation except Figure A1 part e, while results are not satisfactory in pathways and their surroundings. In addition, obvious distortions are apparent in the shadowed areas. Detailed investigation on Figure A3 (e) reveals that there is an obvious distortion in snowy parts of the frame almost in all methods except the proposed CIELab, which resulted in a nearly blue color instead of white snow color. However, CIELab method could be able to preserve the texture and keep the small variances in the color when compared to original MS image. Besides, visual interpretation of the CIELab Pan sharpening results (part c in all figures) demonstrated that use of this color space for Pan sharpening could help to distinguish different tree types and vegetation from each other in the absence of NIR band.
The seven quality metrics, which were presented in Table 2, were used for spectral quality assessment of the Pan sharpening results. Numerical results from these metrics for the rural frames (F1, F5, F6 and F10) are presented in Table A1. The metric values were calculated band by band, and the average values of three bands were used for the accuracy assessment procedure. Numerical results of ERGAS, RASE, RMSE and SAM metrics indicated that the proposed CIELab method produced better results than the remaining methods and was followed by the GS method for most of the metrics. Metric-based results were in line with the visual interpretation. The CIELab method also provided the highest accuracies according to the QAVG and PSNR metrics. In addition, the proposed method provides the value 1 for the SSIM metric, which is the best possible value. Moreover, the IHS method provides worst results for all quality indexes, with respect to Table A1. Lastly, Ehlers, GIHS and HCS methods provide lower accuracies in some cases. This unstable manner of these methods across different scenes is another problem that should be considered.
To assess the spatial quality of Pan sharpened images, the CC, the Zhou index, Sobel RMSE and spatial ERGAS indexes that are presented in Section 3.3, were calculated by comparing the Pan image and the intensity component of the Pan sharpened images. Numerical results of these metrics are presented in Table A2. According to comparative evaluation, the Pan sharpened image from the proposed CIELab method provided the highest spatial CC and Zhou values and lowest SRMSE and SP ERGAS values. These results indicate that the proposed method has the best spatial performance among all methods tested. Ehlers, HIS and HCS methods provided the lowest spatial performances according to the values presented in Table A2.
As a result, the proposed CIELab method provided the best performance for the rural scenes based on the visual interpretation and spectral and spatial quality metrics results. Figure A5 through Figure A7 belongs to the Pan sharpening results of the urban test sites (frames F2, F4 and F9 from D1, D2 and D3 datasets respectively). Each figure belongs to a representative part from the whole image focusing on the buildings and roads, and presents visual comparison between different Pan sharpening techniques on the differently sized and oriented buildings and roads in the urban areas.

Experimental Results from Urban Areas
Visual comparison results of urban areas revealed that all the Pan sharpened images inherited the high spatial information from the Pan image, and likewise, the results of rural areas. Roads and buildings could be better identified in all Pan sharpened images compared to original MS image. As for spectral information, Figure A5 c,e,h, informed us that CIELab, GS and NNDiffuse methods preserved the spectral characteristics and color information in urban areas. In particular, the color information from the buildings with brick rooves are similar to the original MS image. Visual comparison of Figure A5 part a with parts d, f, g and i illustrated that of IHS, HCS, GIHS and Ehlers methods are not able to preserve the original spectral characteristics of buildings as well as the other three approaches did. In particular, Ehlers, HCS and IHS methods provided blurred and smoggy results with faded and paled colors. Parts g and I from Figures A6 and A7 support that HIS and Ehlers methods provide worst visual results among all methods tested. Part e in Figures A6 and A7 reveals the weakest side of GS method; that is, the poor performance in the Pan sharpening of white tones. White colors tend to seem blueish in results of this method. It is obvious from part f in Figures A6 and A7 that the HCS method provided the most blurred result. GIHS and NNDiffuse methods have acceptable results in comparison with the results from other methods (except CIELab method). Detailed investigation of Figures A6 and A7, parts d and h, prove that NNDiffuse method produces distortion in the shadowed areas and GIHS method has poor performance in vegetated areas and trees. Visual interpretation of Figure A5(c), Figure A6(c) and Figure A7(c) reveal the fact that the proposed CIELab method protected spectral properties of original MS image more than the other methods.
Numerical results of spectral quality assessment of Pan sharpened images belonging to the urban test sites (F2, F4 and F9) are presented in Table A3. Metric values demonstrated that the CIELab method provided the most promising results among all Pan sharpening methods used in this research. This method presented the lowest values for the ERGAS, RASE, RMSE and SAM metrics and highest values for QAVG, PSNR and SSIM metrics (again, the highest possible value obtained for SSIM). HCS and IHS methods provided the worst results for most of the metrics. Once again, the second performance rank for spectral quality was obtained by GS method in most of the metrics. Table A4 presents the spatial quality metrics results that were calculated from Pan image and the intensity component of Pan sharpened images for the urban test sites. Similar to the rural test sites, the proposed CIELab method provided the highest CC and Zhou values alongside of lowest SRMSE and SP ERGAS values for urban images, which demonstrated the high spatial quality. In particular, there is great gap between the numeric results of SRMSE and SP ERGAS indexes presented with CIELab method and other methods. Ehlers, IHS and HCS methods acted as the worst methods in the case of spatial indexes, which is consistent with the visual results. Consequently, the proposed CIELab method provided the best performance for the urban test sites (frames F2, F4 and F9) as well, based on the visual interpretation and spectral/spatial quality metrics. Similar to the urban and rural areas, visual comparison of original MS and Pan sharpened images of this category revealed that all the Pan sharpened images produced higher spatial information than original MS image and benefited from Pan image detail level. However, the visual performance of suburban areas was variable, unlike the urban and rural areas. Results from frame F3 ( Figure A8) show that all methods had acceptable performance except Ehlers and HIS. Nevertheless, small amount of distortion in vegetation and shadowed areas is apparent in the results of NNDiffuse and HCS methods. Figure A9 illustrates the effectiveness of the CIELab method in Pan sharpening process. There is an obvious color distortion in all methods except proposed Lab method's result. Parts e, g and i demonstrate similar distortion in the results of Ehlers, GS and HIS methods with a green dominant color distortion, while other three methods, which are presented in parts d, f and h, have purple dominant color distortion. These color distortions are apparent for all surface types including roads, rooves and other objects. The CIELab was the only method that provided acceptable performance for this test frame. Ehlers and IHS methods could not provide good performance for the last test frame, as is observable in parts c and g of the Figure A10. Parts d and f prove that the results of the GIHS and HCS are blurred and not acceptable. Green and white tones distortion is obvious in the result of GS (Part e from Figure A10). The CIELab method illustrates the best performance again in this frame. Regardless of the distortion in shadowed areas, NNDiffuse provided most similar results to original MS image after the CIELab results.

Experimental Results from Suburban Areas
Numerical results of spectral quality assessment of Pan sharpened images belonging to suburban areas (frames F3, F7 and F8) are presented in Table A5. Once again, CIELab method provides the highest values for the QAVG, PSNR and SSIM metrics, while it achieves the lowest values for the ERGAS, RASE, RMSE and SAM metrics. The GS method has the second place again, similar to the urban and rural test sites by achieving better numeric values for most of the metrics. Similar to previous test site results, the Ehlers and IHS methods have the worst performance between all tested methods.
Spatial quality assessment of Pan sharpening results for the frames F3, F7 and F8 are presented in Table A6. The proposed CIELab method illustrated an unrivaled performance in the case of the spatial quality metrics. For the CC and Zhou index, the highest correlation values that indicate high spatial quality are provided by the CIELab method. Moreover, the lowest SRMSE and SP ERGAS values achieved by the proposed method are another proof of high spatial quality of this method. Ehlers, HCS and IHS methods present the lowest CC and Zhou values with the highest SRMSE and SP ERGAS values, which indicate the poor spatial performances of these methods, like their spectral performances.

Thematic Accuracy Evaluation with Spectral Index
The performance of the conventional and proposed Pan sharpening methods were evaluated with several quality indices and visual interpretation in this research. However, the effects of Pan sharpening on the information extraction, process such as image classification, segmentation and index-based analysis is another important concern that requires conservation of spectral properties of the MS image after Pan sharpening. One of the indirect methods frequently used to evaluate the abovementioned situation is to apply the spectral index on the MS and Pan sharpened images, and investigate their consistency. As only visible bands of the images were used in this research, the visible atmospherically resistant index (VARI) proposed by Gitelson et al., 2002 [54] was used for the evaluation as it uses the all visible bands for calculation.
The test site F3 was selected for this evaluation, as it is one of the most challenging sites in the dataset due to complex and heterogeneous land cover characteristics. The VARI index was applied on both the MS and CIELab Pan sharpened images, and a binary classification was performed with the use of the same threshold to map the manmade and natural lands in the region. According to results presented in Figure 4, same level of information extraction could be achieved with CIELab Pan sharpened image and it even provided better thematic representation by providing better geometric representations of the objects and less of the salt and pepper effect observed in the vegetated areas located in the north and south parts of the image.

Overall Comments
When the numerical values from the seven spectral metrics for three different test sites (Tables A1,  A3 and A5) were evaluated, the proposed CIELab had a consistent behavior for different metrics and for different land categories, and ranked as the first among all methods. On the other hand, for the other Pan sharpening methods, different metrics provided various accuracies and did not show a consistent manner for different metrics and even, for different images, the same metric. As an example, the GS method generally had the second ranking for the ERGAS metric for different test sites. However, in the case of the RASE metric, the GS method had second ranking just for some images, while the GIHS method took the second rank for the remaining images. This phenomenon is similar for the worse results; there is no one method that can be mentioned as the worst for all test sites and all metrics. All facts about the consistent manner of the proposed method can also be asserted for the spatial metrics. The CIELab method had the best spatial performance for the all ten test sites, while second through seventh rankings were variable across different metrics and different images. As an outcome, it is evident that the proposed method presents the best results considering spectral and spatial quality metrics and visual interpretation for ten different sites having different landscape characteristics. Moreover, it provided efficient spectral conservation performance according to comparative evaluation performed with binary classification of VARI index. An overall ranking is provided in Table 5, according to expert judgement by considering the quantitative, metric-based results and visual interpretation results together for different landscapes across spectral and spatial domains. Lastly, in order to check the consistency of spectral quality indices through each band of the images, these indices were calculated band-by-band for the test site F3 (Table A7). According to this evaluation, the averaged values for each index are in accordance with the band-based calculations and indices, providing consistent characteristics across image bands in most of the cases.
The test site F3 was selected for this evaluation, as it is one of the most challenging sites in the dataset due to complex and heterogeneous land cover characteristics. The VARI index was applied on both the MS and CIELab Pan sharpened images, and a binary classification was performed with the use of the same threshold to map the manmade and natural lands in the region. According to results presented in Figure 4, same level of information extraction could be achieved with CIELab Pan sharpened image and it even provided better thematic representation by providing better geometric representations of the objects and less of the salt and pepper effect observed in the vegetated areas located in the north and south parts of the image.

Overall Comments
When the numerical values from the seven spectral metrics for three different test sites (Table  A1, A3 and A5) were evaluated, the proposed CIELab had a consistent behavior for different metrics

Conclusions
This research proposed an effective, component substitution-based image Pan sharpening method that uses CIELab color space for Pan sharpening of the VHR Pléiades satellite images. Ten test sites with different landscape characteristics were selected to evaluate the performance of the proposed method in conjunction with six common Pan sharpening algorithms; namely, GS, HCS, IHS, EHLERS, NNDiffuse and GIHS. The comparative evaluation results from Pléiades VHR images supports that the proposed CS algorithm is powerful and ensures better performance compared to the other Pan sharpening methods according to the spectral and spatial accuracy assessment procedures and the visual interpretation. In addition, results indicated that proposed method provided comparatively consistent results, while the performance of other methods varyied with respect to land surface characteristics of the region. As an example for RMSE metric, the best values among the all ten sites were obtained for forest and vegetated areas. Pan sharpening in urban areas resulted in coarser metric values, which illustrate the impact of different land characteristics on the performance of Pan sharpening algorithms. Characteristics of unique CIELab color space, led to producing similar brightness characteristics in Pan sharpened images compared to original MS image. The nature of L* component of MS image helps to preserve spectral and spatial information of original MS and Pan images, respectively. Further improvement of the CIELab-based method could be the implementation of this approach for Pan sharpening of satellite images with more than three bands. In addition, further studies are planned to evaluate the performance of CIELab in fusions of satellite images from different sources. Lastly, other accuracy assessment approaches, such as comparisons of classification and segmentation results of Pan sharpened images, could also help future investigations. Acknowledgments: Great appreciation to ITU CSCRS for providing VHR "Pléiades" images. The authors would also like to acknowledge the many useful contributions of Elif Sertel.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Appl. Sci. 2019, 9,5234 15 of 32                   List A1. Definition of Terms in Table 2 RMSE: MN is the image size, PS(i, j) and MS(i, j) represent pixel digital number (DN) at (i, j) 'th position of Pan-sharpened and MS image. ERGAS: dh dl represents the ratio between the pixel size of high resolution and low resolution images; e.g., 1 4 for Pléiades data, and n number of bands. The RMSE represents root mean square error of band i. SAM: The spectral vector V = {V 1 , V 2 , . . . , V n } stands for reference MS image pixels and V = {V 1 ,V 2 , . . . ,V n } stands for Pan-sharpened image pixels rep reference and both have L components. RASE: The µ represnts the mean of b th band; b is the number of bands and RMSE represents root mean square error. PSNR: The L represents the number of gray levels in the image; MN is the image size, I r (i, j) is pixel value of reference image and I p (i, j) is the pixel value of Pan-sharpened image. A higher PSNR value indicates more similarity between the reference MS and Pan-sharpened images. QAVG: The x and y are the means of reference and Pan-sharpened images, respectively; σ xy is the covariance and σ 2 x and σ 2 y are variances. As QI can only be applied to one band, the average value of three or more bands (QAVG) is used for calculating a global spectral quality index for multi-bands images. QI values range between −1 and 1. A higher value indicates more similarity between reference and Pan-sharpened image. SSIM: The µ stands for mean, σ stands for standard deviation; I r and I p represent reference and Pan-sharpened image respectively. The C1 and C2 are two necessary constants to avoid the index from a division by zero. These constants depend on the dynamic range of the pixel values. A higher value of the measured index shows the better quality of Pan-sharpened algorithm. CC: C r,f is the cross-correlation between reference and fused images, while C r and C f are the correlation coefficients belonging to reference and fused images respectively. SRMSE: Edge magnitude (M) is calculated via spectral distance of horizontal and vertical (M x and M y ) edge intensities. Sp-ERGAS: dh dl represents the ratio between the pixel size of MS and Pan images, and n is the number of bands. Spatial RMSE is represented as below: where MN is the image size, PS(i, j) and Pan(i, j) represents the pixel digital number (DN) at (i, j) 'th position of Pan-sharpened and Pan image.