1.

Introduction

For decades, various image fusion techniques have been developed to obtain a high-resolution multispectral (HRM) satellite image from a set of sensor data: a high-resolution panchromatic (HRP) image containing only intensity information and several low-resolution multispectral (LRM) images with color information.¹ Among them, the intensity-hue-saturation (IHS) domain method and principal component analysis (PCA) method replace the intensity component of LRM images with that of an HRP image and a PCA-applied HRP image, respectively.^1,2 These methods are relatively simple and easy to implement but known to cause distortion in color information.

Recently, wavelet decomposition was adopted for the satellite image fusion, which is reported to have less color distortion.^3–5 In the conventional wavelet-based methods, the “à trous” algorithm is known to be suitable for obtaining wavelet planes of satellite images,^3,6 where the wavelet plane is computed from the difference between the given image and its filtered image with the $B_3$ cubic spline function. The wavelet-based methods are generally classified into the substitution method [substitute wavelet (SW)] and the additive method [additive wavelet (AW)] according to the merging strategy of wavelet planes.³ In the substitution method, the wavelet planes of LRM images are substituted by the wavelet planes of a panchromatic image, whereas those are added by the wavelet plane of a panchromatic image in the additive method. Since the SW discards high frequency components of the LRM images, it usually fails to attain sufficient high frequency details and sometimes loses the information in the LRM images. In contrast, the AW considers the high frequency components of all the LRM images and the HRP image, and thus possibly introduces excessive high frequency details in the synthesized image. Accordingly, there have been efforts to overcome these drawbacks using some improved wavelet-domain methods utilizing weighted merging.^4,5 Otazu et al.⁴ proposed a method in which the wavelet plane of the HRP image is added to each LRM image in proportion to its color intensity value [AW-luminance proportional (AWLP)], whereas Kim et al.⁵ proposed to add the difference between the wavelet planes of HRP image and each LRM image with or without considering the relative radiometric signature of the LRM images [improved AW (IAW) and IAW proportional (IAWP)].

To develop a more optimized way of adding the wavelet planes, we propose a generalized fusion equation in the form of a weighted composition of wavelet planes, where the weights are determined by the least-squares method. The proposed fusion equation includes wavelet planes of HRP, LRM, and degraded HRP images to control the high frequency injection. In the experiments with IKONOS and QuickBird satellite images, the proposed method is compared with the conventional methods in terms of various objective quality metrics. The results show that the proposed method does not introduce noticeable color distortion and enhances the details better than the conventional methods.

2.

Proposed Wavelet-Based Fusion Method

We propose a wavelet-based image fusion method which extracts wavelet planes by the “à trous” algorithm³ and combines them using a generalized fusion equation. The weights included in the fusion equation are determined by the least-squares method.

2.1.

Generalized Fusion Equation

The generalized fusion equation includes $n$ wavelet planes and is defined as

Eq. (1)

$${\mathrm{HRM}}_i={\mathrm{LRM}}_i+{\alpha }_i\sum _{j=1}^n{\omega }_{\mathrm{HRP},j}+{\beta }_i\sum _{j=1}^n{\omega }_{{\mathrm{LRM}}_i,j}+{\gamma }_i\sum _{j=1}^n{\omega }_{\mathrm{LRP},j},$$

where ${\mathrm{HRM}}_i$ is the $i$’th HRM image that we want to obtain, ${\mathrm{LRM}}_i$ is the $i$’th LRM image, and ${\alpha }_i$, ${\beta }_i$, ${\gamma }_i$ are the weights for the $i$’th multispectral image fusion. In this fusion equation, the low-resolution panchromatic (LRP) image is a spatially degraded image of HRP through the decimation (by four) and interpolation (by four), and ${\omega }_I$ is the wavelet plane of the image $I$, e.g., ${\omega }_{\mathrm{HRP}}$ is the wavelet plane of the HRP image. Each LRM image is interpolated (by four) and then added to the weighted sum of the HRP wavelet planes. The weighted LRM wavelet planes and LRP wavelet planes are added or subtracted to inject high frequency components without introducing excessive high frequency components. The “generalized” equation means that our fusion Eq. (1) includes the conventional wavelet-based fusion methods by setting the weights as shown in Table 1.

Table 1

Various wavelet-based fusion algorithms derived from Eq. (1).

	Substitute wavelet 3	Additive wavelet (AW)3	AW-luminance proportional (AWLP)4	Improved AW (IAW)5	IAW proportional (IAWP)5
${\alpha }_i$	1	1	${\mathrm{\Lambda }}_i$	1	${\mathrm{\Lambda }}_i$
${\beta }_i$	$-1$	0	0	0	0
${\gamma }_i$	0	0	0	$-1$	$-{\mathrm{\Lambda }}_i$
${\mathrm{\Lambda }}_i\mathrm{\triangleq }\frac{{\mathrm{LRM}}_i}{(1/n)\sum _{i=1}^n{\mathrm{LRM}}_i}$

As shown in Fig. 1, the LRM images are enlarged to the size of the HRP, and then histogram matching between the HRP and the intensities of the LRM is performed. For the wavelet decomposition of these images, we adopt the “à trous” algorithm as presented in Refs. 3–6. The “à trous” algorithm is known to be more suitable for image fusion than the Mallat algorithm in terms of artifacts and structure distortion, because the “à trous” algorithm is an undecimated and dyadic algorithm which preserves the structure continuity, whereas the Mallat algorithm is a decimated algorithm which causes the loss of linear continuity.^3,6

Fig. 1

The block diagram for the proposed image fusion procedure.

3.

Experimental Results

Twenty-eight IKONOS and 10 QuickBird images were used to evaluate the performance of the proposed fusion method. The spatial resolutions for the HRP and LRM images are 1 and 4 m, respectively. The HRP size is $512\times 512$ and the LRM size is $128\times 128$. It is difficult to measure the quality of the synthesized HRM images objectively since the original HRM images for the reference do not exist. For the quantitative comparison of the fused image, we used the original LRM images with a 4-m resolution as the reference HRM images. Then, the reference HRM images are compared to the HRM images which are obtained by fusing the degraded HRP and LRM images from the original satellite images to 4 and 16-m resolution, respectively.

Figures 2 and 3 show the HRM images synthesized by the proposed method and the existing algorithms, such as IHS,² substitute wavelet intensity (SWI),⁸ and IAWP.⁵ The reference HRM and input LRM images are also included. In these figures, the input LRM images are upsampled by four and interpolated for the purpose of comparison. In Fig. 2, color distortion is noticeable for the IHS method, whereas it is not noticeable for the wavelet-based methods (SWI, IAWP), including the proposed method. It can also be seen that the details are well restored in the synthesized HRM image by the proposed algorithm when compared with the others (Fig. 3). For the objective evaluation, we compute various visual quality metrics,¹ such as correlation coefficient (CC), root mean squared error (RMSE), mean structural similarity (MSSIM), universal image quality index (UIQI), quality nonrequiring reference (QNR), spectral angle mapper (SAM), ERGAS, and peak signal-to-noise ratio (PSNR). For the metrics CC, UIQI, MSSIM, QNR, and PSNR, a larger value means a better performance, whereas a smaller value implies better performance for RMSE, ERGAS, and SAM. The overall visual quality assessments are summarized in Tables 2 and 3. The CPU times for the MATLAB® implementation on a personal computer (Intel Core i5 CPU 750 @2.67 GHz) are also measured for the assessment of the computational complexity. Approximately 80% of the time for the proposed algorithm is occupied by the extraction of the wavelet planes for the HRP, LRM, and LRP images. IHS, Gram-Schmidt adaptive (GSA), GIHSA, AdapIHS, AdapCS, and MMSE are not based on the wavelet decomposition, whereas SWI, AWLP, IAWP, NAW, and the proposed algorithms are based on the wavelet decomposition. In Tables 2 and 3, the best two results for each assessment are highlighted in bold. The assessments show that the proposed algorithm is included in the best two for all the objective evaluation metrics except for the UIQI and QNR tests in the QuickBird images, where the proposed algorithm takes third place for both cases. The performance of the proposed method is comparable to the MMSE method⁹ for the IKONOS data and the proposed method is slightly better than the MMSE except for UIQI in the QuickBird test. The full-size synthesized images for all the algorithms used in the comparison and MATLAB p-codes for the proposed algorithm are available at http://ispl.snu.ac.kr/~idealgod/image_fusion.

Fig. 2

Results of fusing IKONOS1 images [cropped images ($64\times 64$) from fused images ($512\times 512$)]. (a) Reference HRM, (b) input LRM, (c) IHS, (d) SWI, (e) IAWP, (f) proposed method.

Fig. 3

Results of fusing IKONOS2 images [cropped images ($64\times 64$) from fused images ($512\times 512$)]. (a) Reference HRM, (b) input LRM, (c) IHS, (d) SWI, (e) IAWP, (f) proposed method.

Table 2

The spectral quality assessment of image fusion methods (averaged values for 28 IKONOS images). The best two results for each assessment are highlighted in bold.

Table 3

The spectral quality assessment of image fusion methods (averaged values for 10 QuickBird images). The best two results for each assessment are highlighted in bold.

The scaling factor multiplied by the weights is set to 0.65, which is experimentally chosen and applied to all the tested images. Figure 4 shows the performance variation with respect to the scaling factor for IKONOS images. As the scaling factor increases, more high frequency components are injected, whereas a small scaling factor induces less sharpened results. We choose 0.65 as the optimal scaling factor since a scaling factor of 0.65 yields the lowest accumulated ranking values ($=13$) for all the evaluation metrics (cf. 15 for 0.75, 22 for 0.55).

Fig. 4

Performance variation of the proposed method with the scaling factor.

4.

Conclusion

In this article, we have proposed a wavelet domain image fusion method, which can be considered a generalization of the existing methods. The fusion equation consists of the weighted composition of wavelet planes, which includes the terms for enhancing the details and avoiding excessive high frequency components as well. The weights are computed by the least-squares method in the low resolution, and the scaling factor is exploited to compensate the discrepancy incurred in the lower resolution computation of the weights. Experimental results show that the proposed method reduces color distortion and provides high-resolution details.