## 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.^{1} 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}4.^{–}^{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.^{3} 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.^{4} 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.^{5} 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^{3} 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

## (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},$$## 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 | |
---|---|---|---|---|---|

αi | 1 | 1 | Λi | 1 | Λi |

βi | −1 | 0 | 0 | 0 | 0 |

γi | 0 | 0 | 0 | −1 | −Λi |

Λi≜LRMi(1/n)∑i=1nLRMi |

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}

## 2.2.

### Computation of Weights by Least-Squares Method

The fusion equation (1) can be rewritten with matrices as

## (2)

$${\mathrm{HRM}}_{i}=\left[\begin{array}{cccc}{\mathrm{LRM}}_{i}& \sum _{j=1}^{n}{\omega}_{\mathrm{HRP},j}& \sum _{j=1}^{n}{\omega}_{{\mathrm{LRM}}_{i},j}& \sum _{j=1}^{n}{\omega}_{\mathrm{LRP},j}\end{array}\right]{\left[\begin{array}{cccc}1& {\alpha}_{i}& {\beta}_{i}& {\gamma}_{i}\end{array}\right]}^{T},$$## (3)

$${\mathbf{w}}_{i}={({\mathbf{A}}_{i}^{T}{\mathbf{A}}_{i})}^{-1}{\mathbf{A}}_{i}^{T}{\mathbf{Y}}_{i}.$$In the actual implementation, since the HRM image ${\mathbf{A}}_{i}$ is not available, we perform the estimation in the lower resolution. This approach is motivated by an example-based super-resolution, where the relation between an original image patch and the corresponding high-resolution one is extracted using the relation between the original image patch and its degraded (i.e., blurred and downsampled) one.^{7} This process is explained in Fig. 1, which shows that the LRP and LRM images are degraded by decimation and interpolation to generate LLRP and LLRM images (lower resolution images of LRP and LRM), respectively, and ${\mathbf{A}}_{i}$ is built in this lower resolution. In the computation of weights by the least-squares method in the lower resolution, however, the relation between the HRM images and LRM images is possibly different from that between the LRM images and the LLRM images, which might degrade the quality of the fused multispectral images. To compensate for the discrepancy, we adopt a scaling factor which is multiplied by the weights. That is, the weights ${\alpha}_{i}$, ${\beta}_{i}$, and ${\gamma}_{i}$ are computed in the lower resolution and then multiplied by the scaling factor.

## 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,^{2} substitute wavelet intensity (SWI),^{8} and IAWP.^{5} 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,^{1} 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^{9} 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.

## 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.

PSNR | CC | MSSIM | UIQI | QNR | RMSE | ERGAS | SAM | Time (ms) | |
---|---|---|---|---|---|---|---|---|---|

IHS2 | 24.7 | 0.956 | 0.907 | 0.941 | 0.929 | 15.3 | 4.04 | 4.68 | 19 |

SWI8 | 25.9 | 0.963 | 0.901 | 0.947 | 0.967 | 13.3 | 3.49 | 4.94 | 83 |

AWLP4 | 26.6 | 0.968 | 0.919 | 0.955 | 0.965 | 12.7 | 3.36 | 4.63 | 222 |

IAWP5 | 26.6 | 0.968 | 0.920 | 0.955 | 0.964 | 12.7 | 3.35 | 4.63 | 53 |

GSA10 | 25.4 | 0.958 | 0.909 | 0.944 | 0.965 | 14.0 | 3.63 | 4.94 | 32 |

GIHSA10 | 25.5 | 0.965 | 0.908 | 0.944 | 0.968 | 13.7 | 3.63 | 5.00 | 27 |

AdapIHS11 | 25.8 | 0.966 | 0.905 | 0.943 | 0.963 | 13.6 | 3.60 | 4.92 | 593 |

AdapCS12 | 25.5 | 0.962 | 0.888 | 0.935 | 0.956 | 14.5 | 3.79 | 4.81 | 1226 |

NAW13 | 26.1 | 0.964 | 0.905 | 0.948 | 0.967 | 13.1 | 3.45 | 4.94 | 103 |

MMSE9 | 26.7 | 0.969 | 0.931 | 0.964 | 0.969 | 12.0 | 3.15 | 4.71 | 1970 |

Proposed | 26.7 | 0.968 | 0.927 | 0.961 | 0.977 | 12.1 | 3.16 | 4.60 | 463 |

## 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.

PSNR | CC | MSSIM | UIQI | QNR | RMSE | ERGAS | SAM | Time (ms) | |
---|---|---|---|---|---|---|---|---|---|

IHS | 32.4 | 0.934 | 0.926 | 0.952 | 0.948 | 6.52 | 4.35 | 2.98 | 29 |

SWI | 34.1 | 0.951 | 0.937 | 0.958 | 0.967 | 5.52 | 3.72 | 3.00 | 83 |

AWLP | 34.5 | 0.958 | 0.946 | 0.964 | 0.966 | 5.28 | 3.54 | 2.86 | 225 |

IAWP | 34.5 | 0.958 | 0.947 | 0.964 | 0.966 | 5.27 | 3.54 | 2.86 | 51 |

GSA | 32.6 | 0.934 | 0.923 | 0.947 | 0.967 | 6.51 | 4.34 | 3.12 | 33 |

GIHSA | 30.9 | 0.898 | 0.902 | 0.927 | 0.969 | 7.48 | 5.03 | 3.34 | 27 |

AdapIHS | 34.7 | 0.958 | 0.944 | 0.959 | 0.961 | 5.22 | 3.53 | 2.85 | 583 |

AdapCS | 34.5 | 0.959 | 0.939 | 0.953 | 0.937 | 5.55 | 3.81 | 3.21 | 1183 |

NAW | 34.2 | 0.954 | 0.940 | 0.959 | 0.967 | 5.42 | 3.65 | 2.99 | 103 |

MMSE | 34.5 | 0.957 | 0.948 | 0.958 | 0.984 | 5.38 | 3.49 | 2.94 | 1999 |

Proposed | 34.9 | 0.959 | 0.949 | 0.960 | 0.968 | 5.11 | 3.49 | 2.79 | 456 |

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).

## 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.

## Acknowledgments

This work was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT And Future Planning (2009-0083495).