General contourlet pansharpening method using Bayesian inference

18th European Signal Processing Conference (EUSIPCO-2010) Aalborg, Denmark, August 23-27, 2010 GENERAL CONTOURLET PANSHARPENING METHOD USING BAYESIAN INFERENCE Israa Amro1 , Javier Mateos1 and Miguel Vega2 1 Departamento de Ciencias de la Computación e I.A. Universidad de Granada, Granada, Spain iamro@correo.ugr.es, jmd@decsai.ugr.es 2 Departamento de Lenguajes y Sistemas Informáticos Universidad de Granada, Granada, Spain mvega@ugr.es ABSTRACT Pansharpening is a technique that fuses a low resolution multispectral image and a high resolution panchromatic image, to obtain a multispectral image with the spatial resolution and quality of the panchromatic image while preserving spectral information of the multispectral image. In this paper, we present a new pansharpening method based on contourlet transform and Bayesian inference that generalizes the classical contourlet based algorithms. The experimental results show that the proposed method not only enhances the spatial resolution of the pansharpened image, but also preserves the spectral information of the original multispectral image. 1. INTRODUCTION In optical remote sensing, with physical and technological constraints, some satellite sensors provide images of two classes, panchromatic (PAN) images with high spatial and low spectral resolution, and multispectral (MS) images with high spectral and low spatial resolutions. This trade off leads to the lack of high spectral and high spatial resolutions in a single image. The fusion of high spatial resolution PAN images with high spectral resolution MS images creates such an image that is important for many applications, such as feature detection, change monitoring and land cover classification. In general, pansharpening algorithms improve the spatial resolution of the MS image, while simultaneously retaining its spectral information. In the literature, this task has been addressed from different points of view (see [1] and [2] for a description and a comparison of pansharpening methods). In recent years, multiresolution-based methods, which includes Laplacian pyramid [3], wavelet-based methods [4] and contourlet-based methods [5, 6], are becoming popular. The basic idea of all fusion methods based on multiresolution decomposition is to extract the spatial detail information from the PAN image, not present in low resolution MS, to inject it into the later. Non-subsampled contourlet transform (NSCT) [5] is an usual choice in fusion methods, since it provides a complete shift-invariant and multiscale representation, with a fast implementation. Also NSCT can efficiently overcome the limitations of wavelets dealing with high dimensional signals like images (i.e. capture geometric structures in images much more efficiently, because it offers flexible basis functions at more directions and scales). The building block of NSCT are This work has been supported by the Consejerı́a de Innovación, Ciencia y Empresa of the Junta de Andalucı́a under contract P07-TIC-02698. © EURASIP, 2010 ISSN 2076-1465 Figure 1: Non subsampled contourlet transform the 2-D two-channel non-subsampled filter banks (NSFBs). The NSCT is implemented by two shift-invariant stages; a non-subsampled pyramid (NSP) structure that gives the multiscale property, and a non-subsampled direction filter bank (NSDFB) structure that ensures directionality. Both stages of the NSCT, depicted in Fig.1, are constructed to be invertible, in order to have a overall invertible system. Fusion based on NCST can use several methods for the injection of the image details: substitution, the simplest one, addition, and other more complex mathematical models. Regardless of the fusion method, the images to be fused must have the same resolution so a preprocessing step is needed to upsample the MS to the size of the PAN image. In this paper we propose a Bayesian fusion method that comprises, as particular cases, substitution, addition and some complex mathematical models. This paper is organized as follows. In section 2 the general algorithm for NSCT pansharpening, using different injection methods, is described and the used notation is introduced. Section 3 explains the Bayesian modeling and the inference of the high resolution MS image. Experimental results and comparison with other methods are presented in 4 for synthetic and SPOT5 images and, finally, section 5 concludes the paper. 294 2. GENERAL PANSHARPENING ALGORITHM The contourlet pansharpening method is based on the ability of the NSCT transform for obtaining the high frequencies image details in the different directions. Based on this ability, NSCT-based pabsharpening methods intend to extract the detailed information from the PAN image for injecting them into the MS image. The main drawback of NSCT-based pansharpening is the spectral distortion that it may produce. Pansharpening based on multiresolution decompositions, as NSCT, can be carried out in different ways. Following [7] they can be classified as: A. Substitution Model: It involves completely replacing the MS image details, extracted with the NSCT transform, with those of the PAN image. B. Additive Model: Add the NSCT details information of PAN image directly to the MS image bands, or to the NSCT details information of MS image. C. Mathematical-based Model: It is more complex than the above methods since it applies mathematical model to details information in both PAN and MS images and use the model to weight the information of both images in order to control noise and color bleeding effects. All those models start with the observed low resolution MS image, Y, with B bands, Yi , i = 1, . . . , B, each of size P = M × N pixels, and the PAN image, x, of size p = m × n, with M < m and N < n, that contains reflectance data in a single band that covers a wide area of the spectrum. Based on those observations, they find an estimation of y, the high resolution multispectral (HRMS) image, with B bands, yi , i = 1, . . . , B, each of size p = m × n pixels. The general algorithm for pansharpening based on NSCT, is summarized in Algorithm 1, to obtain an estimation of y, ŷ, from x and Y. Algorithm 1 NSCT pansharpening algorithm of x and {Yi } into {ŷi } 1. Upsample each band of the MS image, Yi , to the same size as PAN, x, and register them obtaining si , i = 1, . . . , B. 2. Apply NSCT decomposition on the PAN image x and registered MS image {si }, L x = xR + ∑ D ∑ xld , (1) l=1 d=1 L D si = sRi + ∑ ∑ sldi , i = 1, . . . , B, (2) l=1 d=1 where we are using the superscript R to denote the residual (low pass filtered version) NSCT coefficients band and the superscript ld to refer to the detail bands, with l = 1, . . . , L, representing the scale and d = 1, . . . , D, representing the directions for each coefficient band. 3. Merge the details of PAN {xld } and MS {sld i } images getting {ŷild }, keeping the residual image unchanged, ŷild = az xld + bz sld i , (3) ŷiR = sRi , (4) where z is the decomposition level and directional. 4. Apply the inverse NSCT to merged MS band coefficients {ŷiR }, {ŷild }, getting {ŷi }, L ŷi = ŷiR + ∑ D ∑ ŷild , i = 1, . . . , B. l=1 d=1 (5) Note that for bz = 0, we get the substitution model, and for a = b = 1 we have the additive one, while using different az and bz values we will get different weighted models proposed in the literature. In this paper we propose to modify the merging strategy in step 3 of Algorithm 1 by using Bayesian inference as a mathematical way to estimate the high coefficients of the HRMS image from PAN and MS images and to reconstruct HRMS residual image. The used parameters are estimated at each level of decomposition and direction for each band, providing a sound way to control the noise, preventing color bleeding and generalizing all previous models. 3. BAYESIAN MODELING AND INFERENCE Since the PAN image contains the details of the high resolution MS image but lacks of its spectral information, and the MS image have the spectral information of the HRMS images, the relationship between the HRMS band coefficients and those of the PAN and MS images is defined in this paper as, sRi = yiR + nR , (6) ld + ns i , (7) + nld x, (8) sld i ld x = yild = yild where nR and nld si are the noise of MS residual and coefficients bands, respectively, that is assumed to be Gaussian with zero mean and known variances (βiR )−1 and (βild )−1 , respectively and nld x is the noise of the coefficients bands at each NSCT decomposition level, l, and direction, d, for PAN image, that is assumed to be Gaussian with zero mean and known variance (γild )−1 . Bayesian methods start with the definition of a prior model where we incorporate the expected characteristics of the original NSCT coefficients. Since the residual band is smoothed version of the original HRM band, for yiR we choose a quadratic prior, the Simultaneously Autoregressive (SAR) prior, given by   1 R R (p−1)/2 R 2 R , (9) exp − αi Qyi p(yi ) ∝ (αi ) 2 where Q denotes the Laplacian operator, and αiR is the model parameter that control the smoothness degree of the MS residual band i. For the coefficient band, we choose the Total Variation (TV) prior, that prefers solutions having smooth areas with sharp edges such as the coefficients of the NSCT, given by n o p(yild ) ∝ (αild ) p/2 exp −αild TV (yild ) , (10) q p (∆hk (yild ))2 + (∆vk (yild ))2 where with TV (yild ) = ∑k=1 ∆hk (yild ) and ∆vk (yild ) represent the horizontal and vertical first order differences at pixel k, respectively, and αild is the model parameter of the MS band i coefficients at level l and direction d. Bayesian methods also need of a degradation model, where fidelity to the observed data is incorporated, expressed as the conditional distribution of the observation given the real data. From the observation model of the MS image in 295 Eqs. (6) and (7), we have the following probability distributions   1 2 , (11) p(sRi |yiR ) ∝ (βiR ) p/2 exp − βiR sRi − yiR 2 apply the variational methodology to approximate the posterior distribution by another one, q(yild ), that minimizes the Kullback-Leibler(KL) divergence [8], defined as ld CKL (q(yild )||p(yild |sld i , x )) and =  1 ld ld p/2 ld ) exp − βild sld p(sld |y ) ∝ (β i − yi i i i 2 2  . (12) The probability distribution of the details of the PAN image, given HRM coefficients, from Eq. (8), is written as   2 1 ld ld ld ld p/2 ld ld . (13) p(x |yi ) ∝ (γi ) exp − γi x − yi 2 Having defined the degradation and prior models, the Bayesian inference is performed based on the posterior distribution of the HRMS given the observations. For the residual band, this posterior distribution is given by p(yiR |sRi ) = p(yiR , sRi )/p(sRi ), (14) where p(yiR , sRi ) = p(yiR )p(sRi |yiR ) with p(yiR ) and p(sRi |yiR ) defined in Eqs. (9) and (11), respectively, and the estimation of the real image y, ŷiR , can be calculated as the maximum a posteriori (MAP), that is,   ŷiR = argmax p(yiR |sRi ) = argmin −2 log p(yiR |sRi ) (15) yiR yiR and so h ŷiR = argmin αiR QyiR yiR 2 + βiR sRi − yiR 2 i . (16) (17) where I p is the identity matrix of size p × p. Having found an estimation of the residual bands, let us move to find an estimation of the rest of the NSCT bands. For those coefficients bands, the posterior distribution is formulated as ld ld ld ld ld ld p(yild |sld i , x ) = p(yi , si , x )/p(si , x ) (18) ld ld ld p(sld i |yi )p(x |yi ), q(yild ) )dyild , ld ld ) p(yi |sld , x i (20) which is always non negative and equal to zero only when ld q(yild ) = p(yild |sld i , x ). Unfortunately, the integral in Eq. (20) cannot be directly evaluated due to the TV prior but we can approximate it by using the Majorization-Minimization approach [9] that converts a non-quadratic problem to a quadratic one by the introduction of a new parameter that also needs to be estimated. Thus, the TV prior in Eq. (10) is majorized by the functional (details can be found in [10]), ld p/2 M(yild ,uld i ) = c.(αi )   h (yld ))2 +(∆v (yld ))2 + uld ( j) (∆ i j i j i , q exp −αild ∑ ld j 2 ui ( j) (21) + p ld where uld i is a p−dimensional vector, ui ∈ (R ) , with comld ponents ui ( j), j = 1, ..., p, and that, as we will show later, has a tight relationship with the image. Substituting Eq. (21) into Eq. (10), we obtain p(yild ) ≥ c.M(yild , uld i ). To find an ld estimation of q(yi ) from Eq. (20) we use this lower bound for p(yild ) and alternatively find an estimate of q(yild ) and minimize the bound (see [10] for the details). So, we obtain that h i k+1 h ld 2 v ld 2 (uld ) ( j) = E (∆ (y )) + (∆ (y )) , (22) ld i j i j i qk (y ) for j = 1, . . . , p, where uk+1 ( j) represents the local spatial i ld activity of yi and they will be high for pixels in the neighborhood with high level of detail, thus preserving the image structure, and low for zones with low spatial activity, thus keeping it smooth. By differentiating qk (yild ) = N (yild |Eqk (yld ) [yild ], covqk (yld ) [yild ]), i i (23) we get the estimation of yild , which is the mean of the distribution qk (yild ), where i h k ld ld Eqk (yld ) [yild ] = covqk (yld ) [yild ] (βild (sld i ) + γ xi ) , i i (24) ld ld ld ld ld ld where p(yild , sld i , x ) = p(yi )p(si , x |yi ), with p(yi ) ld defined in Eq. (10) and, assuming that sld i and x are inld dependent, for a given yi , ld ld p(sld i , x |yi ) = q(yild ) log( i By differentiating the right hand side of Eq. (16) with respect to yiR and setting it equal to zero we obtain the estimation of the residual bands ŷiR as yˆiR = (βiR I p + αiR QT Q)−1 βiR sRi , Z ld ld k covqk (yld ) [yild ] = αild ς (uld i ) + βi I p + γi I p i i−1 , (25) with k h t ld k h v t ld k v ς (uld i ) = (∆ ) W (ui ) (∆ ) + (∆ ) W (ui ) (∆ ), (19) ld ld ld where p(sld i |yi ) and p(x |yi ) have been defined in Eqs. (12) and (13), respectively. In order to perform the estimation from the posterior disld tribution we need to calculate p(sld i , x ) in Eq. (18). Howld ld ever, p(si , x ) cannot be calculated analytically and we will h (26) where ∆h and ∆v represent p × p convolution matrices associated with the first order horizontal and vertical differences, respectively, and   k ld k − 21 W ((uld ) ) = diag (u ) ( j) , (27) i i 296 Measure COR Band b1 b2 b3 b1 b2 b3 b1 b2 b3 - SSIM PSNR ERGAS NSCT [6] 0.91 0.91 0.90 0.79 0.81 0.81 26.75 27.17 27.65 5.76 SR [10] 0.84 0.98 0.62 0.90 0.94 0.85 32.68 35.50 30.15 3.12 Proposed 0.99 0.99 0.98 0.97 0.97 0.96 38.17 39.51 36.93 1.61 (a) Original image (b) Observed PAN image (c) Observed MS image (d) NSCT method in [6] (e) SR method in [10] (f) proposed method. Table 1: Synthetic Image Quantative Results is a p × p diagonal matrix, for j = 1, ..., p. This is a spatial adaptivity matrix since it controls the amount of smoothing at each pixel location depending on the strength of the intensity variation that pixel, as expressed by the horizontal and vertical intensity gradient [10]. Returning back to the step 3 in the pansharpening Algorithm 1, and replacing Eqs. (3) and (4) by Eqs. (17) and (25), respectively, we obtain the residual and coefficient bands for the estimated HRMS image. Moreover, by setting αiR = 0, in Eq. (17), we will have the residual band of HRMS estimated image equal to the MS one, that is, it is equivalent to Eq. (3), and setting αild = βild = 0 in Eq. (25) we obtain the NSCT-based substitution model. Also, by setting αild = 0, γild = βild = 1, we get the additive model, while setting gaγild and βild to different values generates different mathematical models. By setting αild to a non-zero value, however, will help our model to control the noise, while merging the coefficients in the different levels and directions. 4. EXPERIMENTAL RESULTS In this section, the proposed NSCT-based pansharpening method using the Bayesian inference is tested. Experiments on a synthetic color image and a real SPOT5 image are conducted to test the proposed method. The observations of the synthetic multispectral are obtained from the color image, displayed in Fig. 2(a), by convolving it with mask 0.25 × 12×2 to simulate sensor integration, and then downsampling it by a factor of two by discarding every other pixel in each direction and adding zero mean Gaussian noise with variance 16. For the PAN image we used the luminance of the original color image and zero mean Gaussian noise of variance 9 was added. The observed MS image was upsampled to the same size of PAN image by the cubic interpolation and then 3 levels of NSCT decomposition was applied on each upsampled MS band and PAN image with 4 and 8 directional levels, for the first two and the third decomposition levels, respectively. The proposed algorithm was run on the resulting coefficients bands until the critek ld k−1 k2 /k(yld )k−1 k2 < 10−4 was satisfied, rion k(yld i i ) − (yi ) k ld where (yi ) denotes the mean of qk (yld i ), which typically is reached within 2 iterations in the first decomposition levels, and 4 iterations in the other levels. The values of parameters were experimentally chosen to be αi = 0.045, βb = 1/16, i = 1, 2, 3 and γ = 0.9. We are working on the method that it can estimate the parameters automatically. We compared the Figure 2: Results for the synthetic image proposed NSCT using Bayesian inference method with the SR method in [10] and the additive NSCT method [6]. The resulted images corresponding to each method are displayed in Fig. 2(d)-(f). To assess the spatial improvement of the pansharpened images we use the correlation of the high frequency components (COR) [1] which takes values between zero and one (the higher the value the better the quality of the pansharpened image). Spectral fidelity was assessed by means of the the peak signal-to-noise ratio (PSNR), the Structural Similarity Index Measure (SSIM) [11], an index ranging from −1 to +1, with +1 corresponding to exactly equal images, and the ERGAS [12] index, a global criterion for what the lower the value, the higher the quality of the pansharpened image. Table 1 shows the corresponding quantitative results. The proposed method provides better results for each measure. The COR values reflect that all methods are able to incorporate the details of the PAN image into the pansharpened one, although the SR method in [10], see Fig. 2(e), introduced less details in the band 3 (blue) since it contributes less to the PAN image and more into the band 2 (green) since it has the highest contribution, which is reflected as a greenish color near the edges of the image. The NSCT method in [6] incorporates details in all the bands but produces a noisy image, see Fig. 2(d). The proposed method (Fig. 2(f)) is able to incorporate details in all the bands while controls the noise and avoids the color bleeding effect. The spectral fidelity measures show that the proposed method performs better than the competing methods, which is also clear from the image in Fig. 2(f). The PSNR for the proposed method is about 297 5. CONCLUSION (a) Observed MS image (b) Observed PAN image In this paper we propose a new pansharpening method that generalizes the fusion strategy of the panchromatic and multispectral images in contourlet based methods. The proposed fusion algorithm is based on the Bayesian modelling and incorporates a solid way to incorporate the details in the panchromatic into the multispectral image while controlling the noise. Particular cases of the proposed fusion method are substitution, additive and weighted contourlets methods. The proposed pansharpening method has been compared with other methods both in synthetic and real images and its performance has been assesed both numerically and visually. REFERENCES (c) Bicubic Interpolation (d) NSCT method in [6] (e) SR method in [10] (f) proposed method. Figure 3: Results for the SPOT5 image 10dB higher than NSCT method in [6] and from 2dB to almost 6dB higher than for the SR method in [10], with a remarkable high SSIM and low ERGAS values which reflect the high quality of the resulting images. In the second experiment, the method was tested on real SPOT5 dataset, where the MS image covers a region of interest of 256 by 256 pixels with pixel resolution of 10 m, while the PAN image is 512 by 512 pixels with a pixel resolution of 5 m. The MS image consists of four bands from the visible and infrared region corresponding to green (0.500.59 µm), red (0.61-0.68 µm), Near IR (0.78-0.89 µm), Mid IR(1.58-1.75 µm), while the PAN image consists of a single band covering the visible and NIR (0.48-0.71 µm). Figure 3(a) shows a region of the RGB color image representing bands 1 to 3 of the MS image. Its corresponding PAN and bicubic interpolation image are depicted in Fig. 3(b) and (c), respectively. The resulted images, displayed in Figs. 3(c)-(f), reveal similar conclusions to the obtained for the synthetic image, from the visual inspections. The NSCT method in [6] (Fig.3(d)) provides a detailed image but quite noisy, the SR method in [10] provides good details for bands 1 and 2, see Fig.3(e), but not for bands 3 and 4 since the PAN image does not cover those bands. This is why the blue color in Fig.3(e), seems to be vanished. The proposed method in Fig.3(f), provides the best result, preserving the spectral properties of MS image while incorporating the high frequencies from the panchromatic image and controlling the noise in the image. [1] V. Vijayaraj. A quantitative analysis of pansharpened images. Master’s thesis, Mississippi State Univ., 2004. [2] L. Alparone, L. Wald, J. Chanussot, P. Gamba, and L. M. Bruce. Comparison of pansharpening algorithms: Outcome of the 2006 GRS-S data-fusion contest. IEEE Trans. Geosci. Remote Sens., 45(10):3012–3021, 2007. [3] B. Aiazzi, L. Alparone, S. Baronti, I. Pippi, and M. Selva. Generalised laplacian pyramid-based fusion of MS + P image data with spectral distortion minimisation. ISPRS Internat. Archives Photogramm. Remote Sensing, 34:36, 2002. [4] M. González-Audı́cana and X. Otazu. Comparison between Mallat’s and the a’trous discrete wavelet transform based algorithms for the fusion of multispectral and panchromatic images. Int. J. Remote Sens., 26(3):595–614, 2005. [5] A. L. da Cunha, J. Zhou, and M. N. Do. The nonsubsampled contourlet transform: Theory, design, and applications. IEEE Trans. Image Process., 15(10):3089– 3101, 2006. [6] M. Lillo-Saavedra and C. Gonzalo. Multispectral images fusion by a joint multidirectional and multiresolution representation. Int. J. Remote Sens., 28(18):4065 – 4079, 2007. [7] K. Amolins, Y. Zhang, and P. Dare. Wavelet based image fusion techniques an introduction, review and comparison. ISPRS. J. Photogramm., 62:249–263, 2007. [8] S. Kullback and R. A. Leibler. On information and sufficiency. Anls. Math. Stat., 22:79–86, 1951. [9] N. Galatsanos. A majorization-minimization approach to total variation reconstruction of super-resolved images. In 16th European Signal Processing Conf., 2008. [10] M. Vega, J. Mateos, R. Molina, and A.K.. Katsaggelos. Super resolution of multispectral images using TV image models. In 2th Int. Conf. on Knowledge-Based and Intelligent Information & Eng. Sys., pages 408– 415, 2008. [11] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli. Image quality assessment: From error visibility to structural similarity. IEEE Trans. Image Process., 13:600 – 612, 2004. [12] L. Wald, T. Ranchin, and M. Mangolini. Fusion of satellite images of different spatial resolutions: Assessing the quality of resulting images. Photogramm. Eng. Remote Sens., 63:691 – 699, 1997. 298