## 1.

## Introduction

Multiresolution image fusion also known as pan-sharpening aims to include spatial information from a high-resolution image, e.g., panchromatic or synthetic aperture radar (SAR) image, into a low-resolution image, e.g., multispectral or hyperspectral image, while preserving spectral properties of a low-resolution image. A large number of algorithms and methods to solve this problem were introduced during the last two decades, which can be divided into two large groups based on a linear spectral transformation followed by a component substitution (CS)^{1} and a spatial frequency decomposition usually performed by means of high-pass filtering (HPF)^{2}^{,}^{3} or multiresolution analysis.^{4}^{,}^{5}^{,}^{6} Sometimes it is quite difficult to orient between all these methods though some classification attempts were already performed.^{6}7.^{–}^{8} We propose to look at these methods from a signal processing view. This type of analysis allowed us to recognize similarities and differences of various methods quite easily and thus perform a systematic classification of most known multiresolution image fusion approaches and methods.^{9} Moreover, this analysis resulted in a general fusion filtering framework (GFF) for a multiresolution image fusion. Sometimes state-of-the-art methods are quite computer time consuming, thus restricting their usage in praxis. To reduce computation time, a simple and fast variant of the GFF method further called a high-pass filtering method (HPFM) is proposed in this paper. It performs filtering in signal domain and thus avoids time-consuming FFT computations.

In parallel to pan-sharpening methods development, many attempts were undertaken to assess quantitatively their quality usually using measures originating from image processing such as mean square error, cross-correlation (CC), structural similarity image (SSIM) index,^{10} Wald’s protocol,^{11} and finally recently proposed joint measures: product of spectral measure based on SSIM and spatial measure based on CC^{12} and quality with no reference (QNR) measure.^{13} Some comparisons are presented in Refs. 14 and 15. Recently a statistical evaluation of most popular pan-sharpening quality assessment measures was performed in Refs. 16, 17, and 18. To measure two different image properties such as spectral and spatial quality, at least two measures are needed, which makes the task of ranking different fusion methods not easy. Following these results, a new joint quality measure (JQM) based both on spectral and spatial quality measures was proposed. This measure was enhanced by proposing a practical normalization of measures and successfully used for a quantitative fusion quality assessment in this paper.

The paper is organized in the following way. First, the four pan-sharpening methods—GFF, HPFM, CS, and Ehlers fusion method—are introduced. Then, a new JQM based on both spectral and spatial quality measures is described. Finally, experiments with a very-high-resolution satellite optical remote sensing WorldView-2 (WV-2) data are performed followed by discussion and conclusions.

## 2.

## Methods

In this section the following four pan-sharpening methods are described: GFF, HPFM, CS, and Ehlers fusion method.

## 2.1.

### General Fusion Filtering

Here the GFF method is shortly introduced in order to better understand the rationale behind a new image fusion method introduced in Sec. 2.2. For detailed description of the GFF method, see Ref. 9.

Let us denote by ${\mathrm{ms}}_{k}$ a low-resolution image, which can be a multispectral/hyperspectral or any other image, with the number of bands equal to $k=1,\dots ,n\text{\hspace{0.17em}}\text{\hspace{0.17em}}(n\in \{1,2,\dots \})$ and by pan a high-resolution image, e.g., panchromatic band, intensity image of synthetic aperture radar (SAR), or any other image. A lot of existing multiresolution methods or algorithms can be seen as an implementation of a general fusion framework:

• Low-resolution image interpolation ${\mathrm{msi}}_{k}=I({\mathrm{ms}}_{k})$;

• Fusion ${\mathrm{msf}}_{k}=F({\mathrm{msi}}_{k},\mathrm{pan})$;

• Histogram matching ${\mathrm{msf}}_{k}=M({\mathrm{msf}}_{k},{\mathrm{ms}}_{k})$.

Each band is processed independently. Further indices are omitted intentionally for the sake of clarity. First and third steps can be included in the fusion step depending on the method. Usually, $I$ is a bilinear (BIL) or cubic convolution (CUB) interpolation and $F$ is a linear or other function of images. In the following we formulate a GFF fusion method including interpolation and fusion in one step.

In order to preserve spectral properties of a low-resolution image ms, one should add only high-frequency information from a high-resolution image pan, thus preventing mixing of low-frequency information of pan with a low-resolution image. The natural way to do it is in a spectral or Fourier domain (signal processing view).

First, both images are transformed into spectral/Fourier space $\mathrm{MS}=\mathrm{FFT}(\mathrm{ms})$ and $\mathrm{PAN}=\mathrm{FFT}(\mathrm{pan})$. Then, high-frequency components are extracted from PAN (solid line, Fig. 1) and added to zero padded spectrum of MS (dotted line, Fig. 1). The formula is written as

where ZP stands for zero padding, $W(f)=\alpha +(1-\alpha )\xb7\mathrm{cos}(2\pi f/{\mathrm{PBW}}_{\mathrm{LR}})$ is Hamming window with $\alpha =0.54$ (aliasing/ringing avoidance), ${\mathrm{PBW}}_{\mathrm{LR}}$ is processing bandwidth of low-resolution image, and HPF is high-pass filter. Cutoff frequency of HPF allows us to control the amount of details added to a low-resolution image. Equivalently we can rewrite Eq. (1) for a low-pass filter (LPF) asFinally, the inverse Fourier transform delivers a fused image with an enhanced spatial resolution

Thus GFF method performs image fusion in Fourier domain. Due to known equivalence formula in signal and spectral domains and assuming that interpolation is performed in a signal domain, we can rewrite Eqs. (1) to (3) as

## (4)

$$\mathrm{msf}=\mathrm{msi}+\mathrm{pan}*\mathrm{hpf}=\mathrm{msi}+\mathrm{pan}-\mathrm{pan}*\mathrm{lpf}.$$Equation (4) defines a fusion function $F$ introduced above and helps to better understand the relation of the proposed method to already known methods.^{9} Here we have to note that Eq. (4) is not the same as GFF due to different interpolation methods used.

GFF method was first introduced in Ref. 9, where its similarity and differences to a high-pass filtering method (Ref. 2) is discussed.

So the only parameter to be selected is a cutoff frequency parameter ${f}_{\text{cutoff}\_\mathrm{HR}}$ controlling the amount of filtering. In this paper a Gaussian filter is used.

## (5)

$$\mathrm{LPF}=\mathrm{exp}[-\frac{1}{2}{\left(\frac{f}{{f}_{\text{cutoff}\_\mathrm{HR}}}\right)}^{2}].$$Of course, any other filter, e.g., Butterworth, can be used. Our experience showed no significant influence of a filter type on the fusion quality. Thus the GFF method depends only on the filtering parameter ${f}_{\text{cutoff}\_\mathrm{HR}}$. We will see in Sec. 4.4 that the optimal parameter for ${f}_{\text{cutoff}\_\mathrm{HR}}$ is $\sim 0.15$ for WV-2 satellite data, which is well supported by other studies (e.g., Ref. 19) and existing experience. Further studies on more data are planned. Moreover, we have to note that the band-dependent filtering parameter can be easily implemented and can increase the quality of pan-sharpening as already known from other studies.^{20}

## 2.2.

### High-Pass Filtering Method

As seen from Eq. (4), it is possible to avoid Fourier transform in the GFF method in order to reduce computation time. Thus a simple and fast variant of the GFF method, further called HPFM, can be derived as follows:

1. Instead of zero-padding in spectral domain, an interpolation of multispectral image in signal domain can be performed.

2. High-pass filter HPF or low-pass filter LPF with desired properties can be built up in a spectral domain and then transformed into signal domain hpf or lpf, respectively.

3. Equation (4) is applied for image fusion in signal domain using convolution with a designed filter.

4. Finally histogram matching is performed.

Proposed method differs from a method introduced in Ref. 2 in the way the filter is constructed and histogram matching performed. Usually a simple boxcar filter is used in signal domain, thus making a fine selection of amount of filtering difficult. For the HPFM the same filters as for GFF method, e.g., Eq. (5), can be used. Here we have to note that due to addition/subtraction operations in Eq. (4), the HPFM can be interpreted as an additive model. Multiplicative model can be written as

Thus the HPFM method depends on the filtering parameter ${f}_{\text{cutoff}\_\mathrm{HR}}$, model, and interpolation method.

## 2.3.

### Component Substitution

CS method is one of the simplest and, maybe, oldest image fusion methods. Here a short description following the recent enhancement of intensity hue saturation transformation method^{1} is given. Under the assumption that msi and pan are highly correlated, one can calculate intensity or mean as

Now CS fusion under the assumption of Eq. (7) can be written as

Similarly as for the HPFM, this way of fusion can be seen as an additive model. Multiplicative model can be written as

Thus a selection of model and interpolation method is required for this method.

## 2.4.

### Ehlers Fusion

For a detailed description of the Ehlers fusion method, see Refs. 3 and 19. Its relation to the GFF method is discussed in Ref. 9. The Ehlers method depends on three parameters: interpolation method (cubic convolution was used as recommended by the authors), cutoff frequency ${f}_{\text{cutoff}\_\mathrm{HR}}$ for high-pass filtering of panchromatic image, and cutoff frequency ${f}_{\text{cutoff}\_I}$ for low-pass filtering of intensity calculated according to Eq. (7). For Ehlers fusion in this paper again Gaussian filter [Eq. (5)] was used for both types of filtering. Finally, we have to note that our implementation of the method further called EhlersX in Xdibias software environment was used. The filtering for this method is implemented in a spectral domain. Validation of our implementation with the original Ehlers software implementation in MATLAB resulted in comparable results. As for other filtering methods GFF and HPFM, the optimal parameters for ${f}_{\text{cutoff}\_\mathrm{HR}}$ and ${f}_{\text{cutoff}\_I}$ are $\sim 0.15$ for WV-2 data, which coincides well with the experience of the authors of the method.

## 3.

## Data

WV-2 satellite remote sensing data over the city of Munich in south Germany were used in our experiments. For scene details see Table 1.

## Table 1

Scene parameters for WorldView-2 data over Munich city.

Parameter | |
---|---|

Image date | July 12, 2010 |

Image time (local) | 10:30:17 |

Mode | PAN+MS |

Look angle | 5.2 deg left |

Product | L2A |

Resolution PAN (m) | 0.5 |

Resolution MS (m) | 2.0 |

## 4.

## Quality Measures

The quality of pan-sharpening is usually measured by spectral and/or spatial quality measures to cover both attributes of a processing result. Measures calculated for the whole image are called global. Window-based measures are calculated in selected areas or, e.g., using sliding window and can distinguish image parts with a different quality. The latter measures are outside the scope of this paper.

## 4.1.

### Spectral Quality

Many spectral quality measures already have been proposed in the literature, e.g., Refs. 5 and 15. Recent comparison^{17}^{,}^{18} showed that the correlation (CORR) between original spectral bands and corresponding low-pass filtered and subsampled pan-sharpened bands is one of the best. It allows us to measure a spectral quality or preservation of a pan-sharpening method for individual bands or by averaging for all bands

## (10)

$$\mathrm{CORR}=\frac{1}{n}\sum _{k}\mathrm{CC}({\mathrm{ms}}_{k},({\mathrm{fms}}_{k}*\mathrm{lpf})\downarrow ).$$It has high values (optimal value is 1) for a good spectral preservation and low values for low spectral characteristics preservation. This measure alone is not able to assess the quality of fusion result, because it is calculated only in reduced image resolution/scale.

## 4.2.

### Spatial Quality

The same investigation^{17}^{,}^{18} showed a preference of SSIM between original panchromatic band and pan-sharpened bands for a spatial quality assessment.

It exhibits high values (optimal value 1) for high spatial quality and low values for low spatial quality. Here we have to note that due to different width of spectral spectra of multispectral and panchromatic data, the CC as proposed in Ref. 12 may not be sufficient because of possible mean and standard deviation differences. SSIM allows us to account for such differences much better.

## 4.3.

### Joint Quality Measure

In an ideal case, pan-sharpening method should exhibit both high spectral and spatial quality measure values. But it is not possible practically, because, e.g., for GFF method (also valid for other filtering methods) different parameters (amount of filtering) lead to different image qualities. Thus a larger high-pass filtering parameter value will lead to a higher spectral quality at the same time reducing spatial quality and vice versa [see Fig 2(b)]. None of the known separate quality measures can fulfill this requirement as a sole measure.^{14} Thus, a JQM could be helpful to achieve optimal parameter selection or best trade-off between spectral and spatial quality or find the best method for a particular application.

One could think of a simple average or product^{12} of both measures: one for spectral and another one for spatial (in this particular case SSIM for spectral quality and CC for spatial quality). In such a way derived quality measure can be easily biased due to different range values of the separate measures. Moreover, CC for spatial quality can be insufficient for data exhibiting different spectral properties as already stated in Sec. 4.2.

In this paper we propose a new JQM that is based on CORR for spectral quality^{10} and SSIM for spatial quality,^{11} measures resulting from discussion in Secs. 4.1 and 4.2.

Due to different ranges of two measures (SSIM is usually lower than CORR as can be seen in Fig 2(b)] we propose here to normalize one of the measures before averaging (producing a joint measure) using a linear scaling transform

## (12)

$${\mathrm{SSIM}}_{\text{norm}}=\frac{\mathrm{SSIM}-{\mathrm{SSIM}}_{\mathrm{min}}}{{\mathrm{SSIM}}_{\mathrm{max}}-{\mathrm{SSIM}}_{\mathrm{min}}}\xb7({\mathrm{CORR}}_{\mathrm{max}}-{\mathrm{CORR}}_{\mathrm{min}})+{\mathrm{CORR}}_{\mathrm{min}},$$Now the averaging of corresponding spectral and normalized spatial measures

delivers a much more meaningful JQM, which is more suitable for parameter selection or comparison of different methods.## 4.4.

### How to Normalize Quality Measures?

For the normalization as proposed in Eq. (12), we need to define four parameters: minimum and maximum values of two variables CORR and SSIM. This can be derived in different ways, e.g., from existing experience (subjective) or experiments with a lot of different methods exhibiting different fusion quality, and thus covering the whole range of spectral and spatial qualities. We follow the latter approach with the following optimizations. Due to diversity of content in an image (remote sensing scene) we propose a data-driven estimation of extreme values. To reduce computation time, we propose to use a single method: a filtering-based method, e.g., HPFM, with two different parameters producing best and worst qualities, e.g., small parameter value (0.05) for high spatial and at the same time low spectral quality and large parameter value (0.7) for low spatial and at the same time high spectral quality [see Fig. 2(b)]. Thus only two runs of a very fast HPFM method deliver the required four parameters. JQM around 0.15 suggests an optimal parameter value for HPFM [Fig. 2(a)].

## 4.5.

### How to Use JQM for Comparison of Different Methods?

JQM introduced in Eq. (13) requires knowledge of normalized ${\mathrm{SSIM}}_{\mathrm{norm}}$, which is not available for all methods. As shown in the previous section, it is even not necessary because these extreme values can be derived using a reference method, e.g., HPFM. Then we can rewrite Eq. (13) using Eq. (12)

where## (15)

$$A=\frac{{\mathrm{CORR}}_{\mathrm{max}}-{\mathrm{CORR}}_{\mathrm{min}}}{{\mathrm{SSIM}}_{\mathrm{max}}-{\mathrm{SSIM}}_{\mathrm{min}}},$$## (16)

$$B=-\frac{{\mathrm{SSIM}}_{\mathrm{min}}\xb7({\mathrm{CORR}}_{\mathrm{max}}-{\mathrm{CORR}}_{\mathrm{min}})}{{\mathrm{SSIM}}_{\mathrm{max}}-{\mathrm{SSIM}}_{\mathrm{min}}}+{\mathrm{CORR}}_{\mathrm{min}}.$$Applying extreme value ranges as proposed above (Sec. 4.4) results in the following $A$ and $B$ values (see Table 2), which are data dependent but are very easy and fast to derive. Thus proposed JQM [Eq. (14)] can be used for estimating quality of any pan-sharpening method on these data.

## Table 2

Estimated values of A and B according to Eqs. (15) and (16) from the proposed extreme range values used for calculation of JQM.

Parameter | Value |
---|---|

A | 0.6786 |

B | 0.42 |

CORRmin | 0.9508 |

CORRmax | 1.0 |

SSIMmin | 0.7822 |

SSIMmax | 0.8547 |

Additional margin of 10% added to these extreme values (avoiding out-of-range values) ensures the appropriate comparison or quality assessment of other pan-sharpening methods. Thus the proposed normalization of quality measures is data dependent, which means it should be performed individually for each image/scene. But this is not a really great drawback because the normalization should be estimated only once per scene and a fast fusion method such as HPFM with different parameter settings can be used.

Here we have to note that there exists one more JQM called QNR (Ref. 13), which is not included in this study, but is one of the topics of the next paper.

## 5.

## Experiments

In this section we shall compare different pan-sharpening methods using JQM as proposed in the previous section. Two more known fusion methods were added to the comparison: Amélioration de la résolution spatiale par injection de structures (ARSIS) (Ref. 21) variant, À Trous wavelet transform model for wavelet fusion^{4} (implementation of Ref. 22) and Gram-Schmidt (GS) spectral sharpening implemented in ENVI Software with an averaging method for low-resolution file calculation and bilinear resampling. Results of quality assessment for various methods are presented in the next two sections by using values of $A$ and $B$ from Table 2. Here we have to note that for the quality assessment only bands spectrally overlapping with panchromatic band are used to assure physically justified results.^{8} Thus, the following three bands were excluded from further analysis: coastal, NIR1, and NIR2.

## 5.1.

### Comparison of Interpolation Methods

In many papers, in addition to the assessment of fusion methods, quality measures for only interpolation (resampling of low-resolution image to high-resolution image) are also referred. Comparison of the four most popular (fast enough for operational applications) interpolation methods is presented in Table 3 and visualized in Fig. 3.

## Table 3

Quality measures of various interpolation methods for WorldView-2 Munich data.

Interpolation method | CORR | SSIM | JQM |
---|---|---|---|

Nearest neighbor | 0.9702 | 0.7860 | 0.9618 |

Zero padding | 0.9781 | 0.7542 | 0.9550 |

Bilinear | 0.9948 | 0.7659 | 0.9673 |

Cubic convolution | 0.9934 | 0.7420 | 0.9585 |

BIL and CUB interpolation methods exhibit the best spectral quality (CORR), which is supported by existing experience. More surprising are results of spatial quality (SSIM). Here the best is nearest neighbor followed by BIL. JQM suggests BIL as interpolation method for pan-sharpening, which agrees well with existing experience. Here we have to note that SSIM is designed for quality assessment of fusion methods whereas interpolation method is not a fusion method and contains no information from panchromatic band. Thus the results of SSIM and JQM for interpolation methods should be treated cautiously and cannot be compared directly with the results of the next section.

## 5.2.

### Comparison of Pan-Sharpening Methods

Comparison of quality and speed of six pan-sharpening methods (some of them with different parameters) is presented in Table 4 and visualized in Fig. 4.

## Table 4

Quality measures and computation time of various pan-sharpening methods for WorldView-2 Munich data on Intel Core 2 Quad CPU Q9450 at 2.66 GHz.

Method | Model | Interpolation method | fcutoff_HR | fcutoff_I | CORR | SSIM | JQM | Time (s) |
---|---|---|---|---|---|---|---|---|

1 GFF | — | Zero padding | 0.15 | — | 0.9782 | 0.8362 | 0.9828 | 85.1 |

2 HPFM | Additive | Bilinear | 0.15 | — | 0.9866 | 0.8337 | 0.9862 | 19.8 |

3 HPFM | Additive | Cubic convolution | 0.15 | — | 0.9873 | 0.8318 | 0.9859 | 22.1 |

4 HPFM | Multiplicative | Bilinear | 0.15 | — | 0.9872 | 0.8359 | 0.9872 | 20.3 |

5 HPFM | Multiplicative | Cubic convolution | 0.15 | — | 0.9878 | 0.8346 | 0.9871 | 22.7 |

6 HPFM | Additive | Bilinear | 0.05 | — | 0.9608 | 0.8447 | 0.9770 | 19.8 |

7 HPFM | Additive | Bilinear | 0.7 | — | 0.9956 | 0.7922 | 0.9766 | 19.8 |

8 CS | Additive | Bilinear | — | — | 0.9406 | 0.8207 | 0.9588 | 19.3 |

9 CS | Multiplicative | Bilinear | — | — | 0.9358 | 0.8310 | 0.9598 | 19.7 |

10 EhlersX | — | Cubic convolution | 0.15 | 0.15 | 0.9450 | 0.8491 | 0.9706 | 160.4 |

11 ARSIS | — | — | — | — | 0.9501 | 0.8663 | 0.9790 | 134.8 |

12 GS | — | Bilinear | — | — | 0.9453 | 0.8192 | 0.9606 | 18.5 |

We see that HPFM method [all four parameter settings (methods 2 to 5)] outperforms its parentage GFF method (method 1). This is because zero padding interpolation smoothens multispectral image to a much greater extent than, e.g., bilinear interpolation (compare CORR values in Table 3). Multiplicative model is better than additive model, whereas the difference between BIL and CUB interpolation methods is negligible. Optimal filtering parameter 0.15 [see Fig. 2(a)] was used for all filtering methods in this paper. Existing experience supports this value even for other data and sensors. Following two parameters 0.05 (method 6) and 0.7 (method 7) of HPFM were used to derive normalizing constants $A$ and $B$, because the results exhibit extreme spectral and spatial qualities as seen in Fig. 4(b). Their JQM is comparable with that of ARSIS method known for its high spatial quality. The lowest JQM as expected exhibit both CS methods and GS sharpening. Ehlers fusion finds its place somewhere between this group and ARSIS. These observations are fully supported by existing experience and visual analysis (Fig. 5).

Additionally computation time of the methods is presented in Table 4 for multispectral image size $1024\times 1024$, panchromatic image size $4096\times 4096$, and for all eight bands of WV-2 Munich data. We see that the proposed HPFM pan-sharpening method is more than four times faster than the parentage GFF fusion method. Moreover, the speed of the proposed method HPFM is comparable with that of classical CS and GS methods. Ehlers and ARSIS fusion methods are about two times slower than GFF and thus are less suitable for operational applications.

## 6.

## Conclusions

A simplified version of a GFF method—a fast, simple, and good HPFM—is introduced with a potential for operational remote sensing applications. It performs filtering in signal domain, thus avoiding time-consuming FFT computations.

A new JQM based on both spectral and spatial quality measures (carefully selected from previous studies) is used to assess the fusion quality of six pan-sharpening methods—GFF, HPFM (with different parameter settings), CS, GS, Ehlers fusion, and ARSIS pan-sharpening methods—on a very high-resolution WV-2 satellite optical remote sensing data. Only spectral bands whose spectral spectrum is overlapping with a spectrum of panchromatic band were used for quality assessment to ensure physically consistent evaluation. The proposed JQM allows a comfortable ranking of different methods using a sole quality measure.

Experiments showed that the HPFM pan-sharpening method exhibits the best fusion quality among several popular methods tested (even better than its parentage method GFF) and at the same time more than four times lower computation time than GFF method. Thus HPFM method, being competitive in speed with known fast methods such as CS and GS but exhibiting a much higher quality, is a good candidate for operational applications. A new quality measure JQM allowed the correct ranking of different pan-sharpening methods, which is consistent with existing experience and visual analysis, thus claiming to be a suitable quality measure for selecting the parameters of a particular fusion method and comparison of different methods.

## Acknowledgments

We would like to thank DigitalGlobe and European Space Imaging for the collection and provision of WorldView-2 scene over the Munich city.

## References

## Biography

**Gintautas Palubinskas** received MS and PhD degrees in mathematics from Vilnius University, Vilnius, Lithuania, in 1981, and Institute of Mathematics and Informatics (IMI), Vilnius, Lithuania, in 1991, respectively. His doctoral dissertation was on spatial image recognition. He was a research scientist at the IMI from 1981 to 1997. From 1993 to 1997, he was a visiting research scientist at German Remote Sensing Data Center, DLR; the Department of Geography, Swansea University, Wales, U.K.; Institute of Navigation, Stuttgart University, Germany; Max-Planck-Institute of Cognitive Neuroscience, Leipzig, Germany. Since 1997, he has been a research scientist at German Remote Sensing Data Center (later Remote Sensing Technology Institute), German Aerospace Center DLR. He is the author/coauthor of about 40 papers published in peer-reviewed journals. Current interests are in image processing, image fusion, classification, change detection, traffic monitoring, data fusion for optical and SAR remote sensing applications.