ADJUSTMENT OF WAVELET DETAILS FOR REDUCING RADIOMETRIC DISTORTIONS IN PAN-SHARPENED PRODUCTS Vladimir BUNTILOV Visiting professor, Mahidol University Faculty of Engineering, Department of Computer Engineering Building #3, RM 6255, 25/25 Phuttamonthon 4 Rd., Nakornprathom, 73170 Thailand E-mail: [email protected]

KEY WORDS: image fusion, pan-sharpening, resolution enhancement, wavelet transform ABSTRACT: The procedure of pan-sharpening aims at the product, which a multispectral instrument would observe if it had as high spatial resolution as the panchromatic one. Usually, the highly resolved spatial information from the panchromatic band is merged with the lower spatially resolved multispectral data. Wavelet-based pan-sharpening methods have shown the ability to preserve radiometric characteristics of the original multispectral data in the sharpened product. However, if the fusing images exhibit substantially different radiometric properties, the resulting image may contain various kinds of defects. This paper presents an approach to minimise distortions in the pan-sharpened image by adjusting the wavelet details of the panchromatic band. The presented scheme was validated using very high resolution optical imagery. 1 INTRODUCTION Due to physical and operational constraints of remote sensing instruments the acquired data might either exhibit high spatial characteristics or maintain a great fidelity of spectral properties of the observed objects. Pan-sharpening techniques aim at providing products, which possess high spatial and spectral accuracy simultaneously (Pohl and van Genderen, 1998). An important factor which must be taken into consideration while developing a fusion method is the difference in relative spectral responses of the panchromatic and the corresponding multispectral bands. If the spectral response of any multispectral sensor has little or no overlap with the spectral response of the panchromatic sensor, then the fused product is highly expected to have different kinds of radiometric and spatial distortions (Thomas et al., 2008, Buntilov and Bretschneider, 2008, Wald et al., 1997). The wavelet transform became a standard tool to ensure the radiometric consistency of the fused product. Nevertheless, as reported by extensive performance evaluations of pan-sharpening methods, no any available algorithm is able to produce the result, which fully satisfies the requirements of both radiometric and spatial fidelity (Thomas et al., 2008, Laporterie-Dejean et al., 2005, Wang et al., 2005). Therefore, the development of a suitable pan-sharpening method is still an open problem. This paper presents a novel algorithm for adjusting the high-frequency wavelet coefficients of the panchromatic band according to the multispectral image. The experiments have demonstrated that fusion methods which use the proposed adjustment technique generate imagery with less spatial artefacts and radiometric distortions than conventional pan-sharpening methods. The remaining part of this paper is organised as follows. Section 2 defines the used notation and discusses theoretical background of pan-sharpening. Section 3 presents the novel method for adjusting the wavelet details. The experiments on performance evaluation of the proposed technique are presented in Section 4. Finally, Section 5 concludes the paper.

2 PROPOSED TECHNIQUE 2.1 Theoretical Background Conventional requirement for most pan-sharpening algorithms is that both images have the same geometric size and are precisely aligned. Hereafter, P AN denotes the panchromatic image, while MS denotes an individual multispectral band, which is upscaled and registered to match P AN. A general wavelet-based fusion process can be described as follows: {W DP AN , W AP AN } = WT N (P AN); {W DM S , W AM S } = WT N (MS) W DF = F (W DP AN , W DM S ); W AF = W AM S F = WT

−1

(1)

(W DF , W AF )

where F denotes the fused image, WT N is the forward wavelet decomposition to the Nth level, W Dx and W Ax are the wavelet details and approximation of the corresponding image, F is the transformation method and WT −1 is the inverse wavelet transform. The main task of developing a wavelet-based pan-sharpening procedure is to construct a proper F . The common assumption is that the high quality spatial information is contained in the highly spatially resolved P AN. Hence, only W DP AN are used to construct W DF . As a result, the initial problem reduced to the task of finding an adjustment function Fadj (W DP AN ) which properly transforms W DP AN according to the radiometric properties of MS. Numerous investigations claim that the defects in the resulting images are resulted from an improper model taken between W DP AN and W DM S (Thomas et al., 2008, Buntilov and Bretschneider, 2008, Wald et al., 1997). 3 PROPOSED ALGORITHM g be an ideal multispectral image which possesses spectral characteristics of the original Let MS N g Consider a MS and spatial resolution of P AN, while W DM are the wavelet details of MS. g S 1 sequence of increasing decomposition levels . At the first level, the difference between W DM g S 1 1 and W DM S would be very significant, since W DM S are mainly resulted from the initial rescaling of the multispectral band, thus not bearing significant spatial information. However, as the decomN N position level increases, W DM and W DM S would become more similar, since the working scale g S becomes coarser. Finally, at a sufficiently deep decomposition level, denoted as the analysis level Na Na Na Na , W DM and W DM S can be considered identical. Let Fadj denotes an adjustment function g S g at the analysis which expresses the relationships between the wavelet details of P AN and MS level of decomposition: Na Na a W DM = Fadj (W DPNAN ) g S

(2)

The proposed method relies on the assumption that local relationships between the wavelet details g remain identical for decomposition levels smaller than Na . In other words, once of P AN and MS the adjustment function has been found for Na , it can be used to calculate the wavelet details of g by transforming the corresponding details of P AN, i.e.: MS Na N W DM = Fadj (W DPNAN ), g S

where N ∈ [1..Na − 1]

(3)

WDMS WDPAN

50 0 -50

0

50

100

150

200

250

300

(a) ratio signal

5 0 -5 0

50

100

150

200

250

300

250

300

(b) multiplicative signal

4000 2000 0 -2000 0

50

100

150

200

(c) adjustment function

2 1 0 0

50

100

150

200

250

300

(d) Na Na Na Na Na Na a Figure 1: (a) 1D profiles of W DPNAN and W DM S , (b) W DM S /W DP AN ratio, (c) W DM S ∗ W DP AN signal, (d) final Fadj calculated using confident areas and clipping

The proposed method uses 2D separable undecimated wavelet decomposition with quadratic spline wavelet function (Mallat and Hwang, 1992). The employed wavelet transform produces N two 2D planes of wavelet details for each level of decomposition, i.e. so-called horizontal HDX N N and vertical V DX wavelet details. Thus, it must be noted that the evaluation of Fadja , the calculations expressed by Equation (3), as well as other procedures of the proposed method are performed N N separately along columns of HDX and rows of V DX , if not explicitly stated otherwise. To calculate the initial approximation of the adjustment function, a multiplicative model is used:

Na Fadj =

Na W DM S Na W DP AN

(4)

Na Na Figure 1(a) depicts the corresponding columns of HDM S and HDP AN , while the ratio signal calNa culated from Equation (4) is shown in Figure 1(b). As can be seen, Fadj exhibit large instabilities Na Na in the areas where the magnitude of HDP AN is small compared to HDM S , e.g. within 200–210 elements region. Na Na Another problem with the initial Fadj occurs in the areas in which the magnitude W DM S is close Na Na to zero. In this case, if the corresponding magnitudes of W DP AN are relatively large, Fadj for this region has negligible values. This suppresses spatial details in the adjustment procedure expressed by Equation 3 and creates unnaturally blurred areas in the fused image. Na Therefore, the model described by Equation (2) is relevant only in the areas where W DM S and Na Na W DM have relatively large magnitudes. It was found that the multiplication signal MU = S Na Na W DM S ∗ W DM S is particularly useful for determining such ”high confidence” areas. As can be seen from Figure 1(c), the magnitude of MU Na drops at the positions where any of the wavelet details has small magnitudes (e.g. elements 190–220, 135–150), while MU Na reveals peaks in

the regions of high confidence (e.g. elements 90–110). Thus, the following procedure is used to Na calculate the refined Fadj : 1. Perform clipping of MU Na by setting its value to 0 for elements in which its magnitude is lower than a specified threshold. Experimentally it was found that a threshold of 0.5 ∗ mean(MU Na ) gives acceptable results. 2. Find the positions of the local extrema of MU Na . 3. For each extremum of MU Na find the neighbourhood within which the magnitude of MU Na is larger than 1/e of the magnitude of the extremum. Na 4. Calculate Fadj for each neighbourhood from Step 3. using Equation (4). Na 5. For regions other than from Step 3. set Fadj = 1. Thus, for regions, in which Equations (4) and (3) are not relevant, a simple select-all fusion rule is applied, i.e. W DPNAN , N ∈ [1..Na − 1] are taken without any modification. This assures that ”non-confident” areas are still sharpened, although, perhaps not in the most optimal way.

Na The procedure described above eliminates most instabilities found in the initial Fadj . However, it was found that the fused product may still exhibit some spatial artefacts due to the inconsistent Na Na has erroneously large or small . The first type of problematic areas is where Fadj values of Fadj values in spite of the undertaken measures. Experimentally it was found that for data described in Na Na | < 0.2 are very likely to be improperly calculated. > 5 and |Fadj Section 4 the values of Fadj Na is negative. In this case, although the procedures The other type of artefacts shows up when Fadj correctly identified that the region contains a spatial object which produces contrast inversion, the model described by Equation (3) fails to produce a satisfactory result. Therefore the areas Na Na Na is set to 1 at the < 0.2 are considered as ”non-confident” and Fadj > 5 and Fadj in which Fadj Na corresponding locations. The final Fadj is depicted in Figure 1(d).

4 EXPERIMENTS To validate the proposed adjustment algorithm, real optical imagery from IKONOS-2 satellite was used. MS image was upscaled and registered to match P AN. The following pan-sharpening methods were employed in performance comparison1: 1. Wavelet-based fusion with full details substitution rule (WT) (Yocky, 1995): Each multispectral band was sharpened individually by combining the approximation image with the wavelet details of P AN. The same wavelet as described in Section 3 was used for analysis. The ratio of spatial resolutions of MS and P AN is four, thus the deepest decomposition level was choosen to be two. 2. Wavelet-based fusion based on the proposed adjustment method (WTA): Each multispectral band was sharpened individually by combining the approximation image with the wavelet details constructed as described in Section 3. The third decomposition level was chosen as the analysis level. 3. PCA-based method (PCA) (Shettigara, 1992): In this method the first principal component P C1 of the entire set of the multispectral bands was replaced by the panchromatic band. 1 The

GNU Octave source code of the experiments can be found at http://sites.google.com/site/buntilov/download

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 2: Examples of pan-sharpening of IKONOS-2 imagery. (a) 200 × 200 pixels subscene of colour combination of original M S bands, (b) corresponding P AN image, (c) WT, (d) WTA, (e) PCA, (f) IHS, (g) PCA-WT, (h) PCA-WTA, (i) IHS-WT, (j) IHS-WTA

4. IHS-based method (IHS) (Shettigara, 1992): In this method the imagery was sharpened by replacing its intensity component I by the panchromatic band. 5. Combination of PCA- and wavelet-based techniques (PCA-WT) (Li et al., 1999): In this method the first principal component P C1 of the decorrelated MS bands was replaced by its sharpened counterpart P C1∗ . To obtain P C1∗ , P C1 was spatially enhanced by P AN using the WT method. 6. Combination of PCA- and the proposed techniques (PCA-WTA): The first principal component P C1 of the decorrelated MS bands was replaced by its sharpened counterpart P C1∗ . To obtain P C1∗ , P C1 was spatially enhanced by P AN using the WTA method. 7. Combination of IHS- and wavelet-based techniques (IHS-WT) (Kumar et al., 2000): In this method the intensity component was sharpened by P AN using the WT method. 8. Combination of IHS- and the proposed techniques (IHS-WTA): In this method the intensity component was sharpened by P AN using the WTA method.

Figure 2 depicts the extracted patches of P AN, colour composite of MS bands as well as fused images. The patches were specifically selected to show the cases when WTA outperforms WT. As can be observed, all wavelet-based methods retain radiometric properties of MS relatively well. Contrary, the results of IHS and PCA techniques have significant colour distortions. Visual analysis shows that the techniques, which use the proposed wavelet adjustment method (WTA, PCA-WTA, IHS-WTA), minimise distortions in the areas where P AN and MS have significant differences in spectral properties. For example, consider a rhombus feature in the bottom part and a square object in the upper part of the scene. The images pan-sharpened by common wavelet-based methods (WT, PCA-WT, IHS-WT) demonstrates a ringing problem near the borders in the form of light edging inside the objects. At the same time, the techniques based on the proposed method substantially reduces this artefact.

5 CONCLUSIONS This paper proposes a novel method for adjusting the wavelet details of the panchromatic band according to the multispectral data. The adjustment procedure ensures the radiometric and maintains spatial consistency of the pan-sharpened product. The proposed method was compared with a number of conventional pan-sharpening algorithms using very highly resolved optical data from IKONOS-2 satellite. The experiments demonstrated that the new method is able to reduce the distortions in the fused product. ACKNOWLEDGEMENTS c The imagery for the experiments was downloaded from Space Imaging LLC website: http://www.geoeye.com/CorpSite/resource/sample_imagery.aspx.

REFERENCES Buntilov, V. and Bretschneider, T., 2008. Investigation on image fusion of remotely sensed images with substantially different spectral properties. Vol. 7109, Number 1, SPIE, p. 710903. Kumar, A., Kartikeyan, B. and Majumdar, K., 2000. Band sharpening of IRS-multispectral imagery by cubic spline wavelets. International Journal of Remote Sensing 21(3), pp. 581–594. Laporterie-Dejean, F., de Boissezon, H., Flouzat, G. and Lefevre-Fonollosa, M.-J., 2005. Thematic and statistical evaluations of five panchromatic/multispectral fusion methods on simulated PLEIADES-HR images. Information Fusion 6(3), pp. 193–212. Li, J., Zhou, Y. and Li, D., 1999. PCA and wavelet transform for fusing panchromatic and multispectral images. Proceedings of the SPIE 3719, pp. 369–377. Mallat, S. and Hwang, W. L., 1992. Singularity detection and processing with wavelets. IEEE Transactions on Information Theory 38, pp. 617–643. Pohl, C. and van Genderen, J. L., 1998. Multisensor image fusion in remote sensing: concepts, methods and applications. International Journal of Remote Sensing 19(5), pp. 823–854. Shettigara, V., 1992. A generalized component substitution technique for spatial enhancement of multispectral images using a higher resolution data set. Photogrammetric Engineering and Remote Sensing 58(5), pp. 561–567. Thomas, C., Ranchin, T., Wald, L. and Chanussot, J., 2008. Synthesis of multispectral images to high spatial resolution: A critical review of fusion methods based on remote sensing physics. Geoscience and Remote Sensing, IEEE Transactions on 46(5), pp. 1301–1312. Wald, L., Ranchin, T. and Mangolini, M., 1997. Fusion of satellite images of different spatial resolutions: assessing the quality of resulting images. Photogrammetric Engineering and Remote Sensing 63(6), pp. 691–699. Wang, Z., Ziou, D., Armenakis, C., Li, D. and Li, Q., 2005. A comparative analysis of image fusion methods. IEEE Transactions on Geoscience and Remote Sensing 43(6), pp. 1391–1402. Yocky, D., 1995. Image merging and data fusion by means of the discrete two-dimensional wavelet transform. Journal of the Optical Society of America A 12(9), pp. 1834–1841.