19th European Signal Processing Conference (EUSIPCO 2011) Barcelona, Spain, August 29 - September 2, 2011 PARAMETER ESTIMATION IN THE GENERAL CONTOURLET PANSHARPENING METHOD USING BAYESIAN INFERENCE Israa Amro1 , Javier Mateos1 and Miguel Vega2 1 2 Departamento de Ciencias de la Computación e I.A., Universidad de Granada, Granada, Spain iamro@correo.ugr.es, jmd@decsai.ugr.es Departamento de Lenguajes y Sistemas Informáticos, Universidad de Granada, Granada, Spain mvega@ugr.es ABSTRACT This paper solves the problem of parameter estimation for the general contourlet pansharpening method using Bayesian inference. In the general contourlet panshapening method, a set of parameters that control the contribution of each band of the multispectral image, the panchromatic image and the prior knowledge on the image need to be set. The proposed method takes into account the relationship between contourlet coefficients to incorporate prior knowledge on the unknown parameters in the form of hyperprior distributions. This method is able to estimate all the unknown parameters together with the high resolution multispectral image in a fully automatic way. 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 Pansharpening is a pixel based technique that fuses a low resolution (LR) multispectral (MS) image and a high resolution (HR) panchromatic (PAN) image to provide a HR MS image with the spatial resolution and quality of PAN image and the spectral resolution of the MS image. In the literature, this task has been addressed from different points of view (see [1] for a complete review). In recent years, multiresolution-based methods, which includes Laplacian pyramid, wavelet-based methods and contourlet-based methods, are becoming popular. Non-subsampled contourlet transform (NSCT) [2] is an usual choice in fusion methods, since it provides a complete shift-invariant and multiscale representation, with a fast implementation. Fusion based on NSCT can use several methods for the injection of the image details. Recently, we proposed a general NSCT pansharpening method using Bayesian inference [3] that comprises, as particular cases of injection methods, substitution, addition and some complex mathematical models, but needed of a big amount of parameters to be set. In this paper, we extend the work in [3] to automatically estimate the unknown parameters, within the hierarchical Bayesian framework. This paper is organized as follows. In section 2, the used notation is introduced and the general algorithm for NSCT pansharpening, using different injection methods, is described and the relationship between the NSCT coefficients is studied. Section 3 explains the Bayesian modeling and section 4 presents the inference of the high resolution MS image and the estimation of the parameters. Experimental results and comparison with other methods are presented in section 5 for synthetic and SPOT5 images and, finally, section 6 concludes the paper. 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, 2011 - ISSN 2076-1465 2. GENERAL PANSHARPENING ALGORITHM AND BAYESIAN FORMULATION Contourlet based pansharpening algorithms start with the observed LR MS image, Y , with B bands, Yb , b = 1, . . . , B, each of size P = M × N pixels, which provides low spatial but high spectral resolution, and the PAN image, x, of size p = m × n, with M < m and N < n, a high spatial resolution image which contains reflectance data in a single band that covers a wide area of the spectrum. Based on those observations, contourlet based pansharpening algorithms find an estimation of y, the HR MS image, with B bands, yb , b = 1, . . . , B, each of size p = m × n pixels. Contourlet pansharpening methods are based on the ability of the NSCT transform for obtaining the high frequencies image details at different scales and different directions to extract the detail 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. In general, the pansharpening algorithm in Algorithm 1 can be used to obtain an estimation of y, ŷ, from x and Y . Algorithm 1 NSCT pansharpening algorithm of x and {Yb } into {ŷb } 1. Upsample each band of the MS image, Yb , to the size of the PAN, x, and register them obtaining sb , b = 1, . . . , B. 2. Apply NSCT decomposition on the PAN image x and registered MS image {sb }, x = xr + D L X X xld , (1) l=1 d=1 sb = srb + D L X X sld b , b = 1, . . . , B, (2) l=1 d=1 where the superscript r denotes the residual (low pass filtered version) NSCT coefficients band and the superscript ld refers to the detail bands, with l = 1, . . . , L, representing the scale and d = 1, . . . , D, representing the direction for each coefficient band. 3. Merge the details of PAN {xld } and MS {sld b } images getting {ŷbld }, keeping the residual image unchanged, ŷbld = ald xld + bld sld (3) b , ŷbr = srb . (4) 4. Apply the inverse NSCT to merge the MS band coefficients {ŷbr }, {ŷbld }, getting {ŷb }, D L X X ŷbld , b = 1, . . . , B. (5) ŷb = ŷbr + l=1 d=1 1130 Figure 3: Marginal distribution of the difference between neighbor coefficients and the distribution of TV prior. Figure 1: NSCT coefficient relationships. Starting from modeling the relation between yb and sb , as s b = yb + n b , (a) (b) Figure 2: Conditional distribution of a finest subband of the first band of synthetic image, conditioned on (a) parent p(X|P X) and (b) neighbor p(X|N X) Note that for bld = 0, we get the substitution model [4], and for ald = bld = 1, ∀l = 1, . . . , L and d = 1, . . . , D we have the additive one [5], while using different ald and bld values we will get different weighted models proposed in the literature. In [3] we proposed to modify the merging strategy in step 3 of Algorithm 1 by using Bayesian inference as a mathematical way to estimate the details coefficients of the HRMS image from those of the PAN and MS images. In this paper we extend the formulation in [3] to deal with the problem of parameter estimation. Using the hierarchical Bayesian paradigm we include information on the unknown parameters in the form of hyperprior distributions and estimate the unknown parameters together with the high resolution multispectral image. The parameters incorporated in our model 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. Before going into the formulation of the problem let us study the characteristics of the NSCT coefficients. These characteristics for the subsampled case have been studied in [6], they found three relationships in contourlet coefficients, which are shown in Figure 1. The reference coefficient has eight neighbors (NX) in the same subband, parent (PX) at the same spatial location in the immediately coarser scale and cousins (CX) at the same scale and spatial location but in different directional subbands. Following [6], we studied the conditional distributions of contourlet coefficients, conditioned on their parents, neighbors and cousins, in order to know the relationship between these coefficients for the NSCT. Figure 2 shows the conditional distributions for the parents and neighbors. The distribution for the cousins is not shown due to space constraints but is similar to the other two. We notice that all of these conditional distributions exhibit a ”bow-tie” shape. Therefore, we conclude that contourlet coefficients of natural images are approximately uncorrelated yet dependent on each other. These dependencies, however are local. This means that we can estimate each band of the NSCT coefficients independently but we should take into account the relation between the bands. We are going to relate the different bands through the value of the parameters. In this paper we are going to relate the HR MS image and the upsampled MS image following the model in [3]. (6) with nb being the capture noise assumed to be Gaussian with zero mean and covariance matrix σb2 I. Applying the NSCT to both sides of Eq. (6) and since the NSCT decomposition is a tight frame, we can write D L X D L X X X r r (7) ybld + nld srb + sld b . b = yb + n b + l=1 d=1 l=1 d=1 Assuming that the noise is separable, that is, it decomposes in the same way as the image does, and therefore low frequencies of the noise affect only to the low frequencies of the MS band, and the high frequencies in each direction affect only to their corresponding high frequencies in each direction, we can write that ld ld sld (8) b = yb + nb , for all l = 1...L, d = 1...D. We are going to approximate the noise of the detail bands, nld b , by a Gaussian distribution with zero mean and covariance matrix (βbld )−1 I, with βbld the inverse of the unknown noise variance of the detail band at level l and direction d of the MS band b. Since the PAN image contains the details of the high resolution MS image but lacks of its spectral information, the relationship between the HRMS band coefficients and the PAN image could be written as xld = ybld + nld x , (9) where nld x is the noise of the coefficients bands at each NSCT decomposition level, l, and direction, d, for PAN image, that it is approximated by a Gaussian distribution with zero mean and covariance matrix (γbld )−1 I. Note that, with this modeling, we have decoupled each one of the bands of the contourlet transform and, since they are uncorrelated, we can do the estimation of each band independently of the other bands. 3. BAYESIAN MODELING When modeling the coefficient bands into the hierarchical Bayesian framework we have two stages. In the first stage, knowledge about the structural form of the noise coefficients and the structural behavior of the HR MS imld ld ld age coefficients is used in forming p(sld b , x |yb , Ωb ) and p(ybld |Ωld ), respectively. These noise and image models b depend on the unknown parameters Ωld that need to be b estimated. In the second stage a hyperprior on the parameters is defined, thus allowing the incorporation of information about these hyperparameters into the process. Thus, according to Bayesian paradigm, we can define the joint distribution on the observation, hyper-parameters, ld ld ld and HR MS coefficient bands as, p(Ωld b , yb , x , s b ) = ld ld ld ld ld ld ld ld ld p(Ωb )p(yb |Ωb )p(sb |yb , Ωb )p(x |yb , Ωb ), and inference is ld ld ld based on p(Ωld b , yb |sb , x ). From the observation model of 1131 the MS and PAN image in Eqs. (8) and (9) respectively, we have the following probability distributions   2 1 ld ld ld ld ld p/2 , (10) s − y β p(sld |y ) ∝ (β ) exp − b b b b b b 2   2 1 p(xld |ybld ) ∝ (γbld )p/2 exp − γbld xld − ybld . (11) 2 Following [3, 7], we choose a prior model based on the Total Variation (TV) for the HR MS image, y, coefficient bands given by n o p(ybld |αbld ) ∝ (αbld )p/2 exp −αbld T V (ybld ) , (12) Pp p h ld 2 (∆i (yb )) + (∆vi (ybld ))2 where with T V (ybld ) = i=1 h ld v ld ∆i (yb ) and ∆i (yb ) represent the horizontal and vertical first order differences at pixel i, respectively, and αbld is the model parameter of the MS band b coefficients at level l and direction d. The idea behind this model is to consider the image as a set of relatively smooth objects or regions separated by strong edges, such as the coefficients of the NSCT. As can be observed in Figure 3, the TV prior almost perfectly fits the distribution of the NSCT coefficients for the details bands. In the second stage of the hierarchical Bayesian framework we define the distribution on the parameters. To model the hyperparameters we use a gamma distributions, p(w|aw , cw ) = Γ(w|aw , cw ), (13) ld ld ld where w > 0, w ∈ Ωld b = (αb , βb , γb ) denotes a hyperparameter, and aw > 0 and cw > 0 are, respectively, the shape and the inverse scale parameters of the distribution. Finally, combining the first and second stage of the problem modeling we have the global distribution, ld ld ld ld ld ld ld ld p(Ωld b , yb , x , sb ) = p(αb )p(βb )p(γb )p(yb |αb ) ld ld ld ld ld p(sld b |yb , βb )p(x |yb , γb ), ld p/2 M (αbld , ybld , uld b ) = c.(αb ) " # p X (∆hi (ybld ))2 +(∆vi (ybld ))2 + (uld ld b )(i) p exp −αb , 2 (uld b )(i) i=1 ld + p where uld b is a p−dimensional vector, ub ∈ (R ) , whose ld components (ub )(i), i = 1, ..., p, are quantities that need to be computed and have an intuitive interpretation related to the unknown images ybld . Comparing Eq. (16) with Eq. (12), we obtain p(ybld , αbld ) ≥ c.M (αbld , ybld , uld b ). This leads to the following lower bound for the joint probability distribution ld ld ld ld ld ld ld ld ld ld p(Ωld b , yb , sb , x ) ≥ c.p(Ωb )M (αb , yb , ub )p(sb |βb , yb ) ld ld ld ld × p(xld |γbld , ybld ) = F (Ωld b , yb , sb , x , ub ). (17) Hence, a sequence can be obtained of ever decreasing upper bounds for the KL divergence as ld ld ld ld ld CKL (q(Ωld b , yb )||p(Ωb , yb |sb , x )) ld ld ld ld ld ld k ≤CKL(q(Ωld b ).q(yb )||F (Ωb , yb , sb , x , (ub ) )). (18) The following algorithm which extends the one presented in [3] to deal with parameter estimation, can be used for ld calculating the approximating posteriors q(Ωld b ) and q(yb ). Algorithm 2 Posterior parameter and image distribuld tions estimation in TV reconstruction using q(Ωld b , yb ) = ld q(Ωld )q(y ) b b Given u1 ∈ (R+ )p and q 1 (Ωld b ), an initial estimate of the distribution q(Ωld b ), For k = 1, 2, ..., until a stopping criterion is met. 1. Find q k (ybld ) = arg min q(ybld ) (14) ld ld ld ld ld ld k CKL(q k(Ωld b )q(yb )||F (Ωb ,yb ,sb ,x ,(ub ) )) (19) k+1 2. Find (uld = arg min b ) ld ld ld ld ld where p(ybld |αbld ), p(sld b |yb , βb )) and p(x |yb , γb ) are given in Eqs. (12), (10), and (11), respectively. uld b k ld ld ld ld ld ld CKL(q k (Ωld b )q (yb )||F (Ωb ,yb ,sb ,x ,ub )) (20) 4. BAYESIAN INFERENCE For our selection of hyperparameters in the previld ous section, the set of all unknowns is (Ωld = b , yb ) ld ld ld ld (αb , βb , γb , yb ). The Bayesian paradigm dictates that inld ld ld ld ld ference on (Ωld b , yb ) should be based on p(Ωb , yb |sb , x ) = ld ld ld ld ld ld ld ld p(Ωb , yb , sb , x )/p(sb , x ). Since p(sb , x ) cannot be calld ld ld culated analytically, then p(Ωld b , yb |sb , x ) can not be found in closed form, so, we apply the variational methodology to approximate the posterior distribution by another distribuld tion, q(Ωld b , yb ), that minimizes the Kullback-Leibler(KL) divergence, defined as ld ld ld ld ld CKL (q(Ωld b , yb )||p(Ωb , yb |sb , x ))   ZZ ld q(Ωld ld ld b , yb ) = q(Ωld , y )log dΩld b b b dyb , ld ld ld p(Ωld b , yb |sb , x ) (16) (15) which is always non negative and equal to zero only when ld ld ld ld ld q(Ωld b , yb ) = p(Ωb , yb |sb , x ). We choose to approximate the posterior distribuld ld ld ld ld tion p(Ωld b , yb |sb , x ) by the distribution q(Ωb , yb ) = ld ld ld ld q(Ωb )q(yb ), where q(yb ) and q(Ωb ) denote distributions on ybld and Ωld b , respectively. Due to the form of the TV prior, the above integral cannot be directly evaluated so, following [3, 7], we approximate it by using the MajorizationMinimization approach [8]. Thus, the TV prior in Eq. (12) is majorized by the functional, 3. Find q k+1 (Ωld b ) = arg min q(Ωld ) b k+1 ld ld ld ld ld k+1 CKL(q(Ωld (yb )||F (Ωld )) (21) b )q b ,yb ,sb ,x ,(ub ) k ld ld k ld Set q(Ωld b ) = limk→∞ q (Ωb ), q(yb ) = limk→∞ q (yb ). In Eq. (19), we obtain for q k (ybld ), the p-dimensional Gaussian distribution q k (ybld ) = N (ybld |E k [ybld ], cov k [ybld ]), with h i ld ld , (22) E k [ybld ] = cov k [ybld ] βbld sld b +γ x h i−1 k ld ld cov k [ybld ] = αbld ς(uld , (23) b ) + βb Ip + γ b Ip k k h v t ld k v = (∆h )t W (uld and ς(uld b ) b ) (∆ ) + (∆ ) W (ub ) (∆ ) k where W ((uld ) ) is an p × p diagonal matrix of the form b  − 1  2 k k (uld W (uld , (24) b ) (i) b ) = diag and ∆h and ∆v represent the p × p convolution matrices associated to the first order horizontal and vertical differences, respectively. When an estimation ŷbld of the band coefficients is needed in step 3 of Algorithm 1, the mean of the distribution q(ybld , E[ybld ], is selected. k+1 To calculate (uld from Eq. (20), we obtain b ) 1132 Measure COR Band b1 b2 b3 b1 b2 b3 b1 b2 b3 - NSCT [5] 0.91 0.91 0.90 0.79 0.81 0.81 26.75 27.17 27.65 5.76 SR [7] 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 5. 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 SSIM conducted to test the proposed method. The observations of the synthetic multispectral are obtained from the color image, displayed in Fig. 4(a), by convolving it with mask PSNR 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 lumiERGAS nance of the original color image and zero mean Gaussian noise of variance 9 was added. The observed MS image was Table 1: Synthetic Image Quantative Results upsampled to the same size of PAN image by bicubic inh i k+1 (uld (i) = Eqk (yld ) (∆hi (ybld ))2 + (∆vi (ybld ))2 , i = 1, · · · , p. terpolation and then 3 levels of NSCT decomposition was b ) b applied on each upsampled MS band and PAN image with 4 and 8 directional levels, for the first two and the third k ld ld k+1 Once we know q (yb ) and (ub ) , the next step is to calcudecomposition levels, respectively. The proposed algorithm late the distributions on the hyperparameters from Eq. (21). was run on the resulting coefficients bands until the critek ld k−1 2 k−1 2 We first note that rion k(yld k /k(yld k < 10−4 was satisfied, i ) − (yi ) i ) ld k k where (yi ) denotes the mean of q (yild ). When all the pak+1 q k+1 (Ωld (αbld )q k+1 (βbld )q k+1 (γbld ) b ) = q rameters are automatically estimated with no prior knowledge about them, that is, λw = 0, w ∈ Ω, Algorithm 2 where q(ω) = Γ(ω|āω , c̄ω ) which produces typically converges within 4 iterations. Since, as previously λαld said, the different level and direction bands are indepen(E k+1 [αbld ])−1 = b + (1 − λαld ) dent though related, in order to provide the method with b ld αb some prior knowledge about the value of the parameters, we p   2  2  have studied the relationships between the parameters on 2X k ld k ld v h × + ∆i E [yb ] ∆i E [yb ] the different levels. For that, we generated images following p i=1 a Gaussian distribution with zero mean and different values  1 for the variance, applied the contourlet decomposition to the h i 2 1 images, and calculated the variance of the noise at each level + tr (cov k [ybld ])−1 ×((∆h )t(∆h )+(∆v )t(∆v )) p and direction. Studying those values we found the relation(25) ship between the inverse of the variance of the coefficients for three level of decompositions, the first with 4 directions λ ld β (E k+1 [βbld ])−1 = b + (1 − λβ ld ) and the second and third with 8 directions, those values were b βbld used for β¯bld and γ¯bld and we set λ1d ω = 0, for all ω ∈ Ω, and
ld ld ld ld ld 2 k ld −1 λ sb − Eqk (y) [yb ] +tr (cov [yb ]) αb = 0, λβb = 0.9, λγb = 0.9, for l > 1. We found that , (26) providing the algorithm with this prior knowledge, from the p relationship between the coefficients, improved the obtained λγ ld images and reduced the calculation time. We compared the k+1 ld −1 b + (1 − λγ ld ) (E [γb ]) = proposed NSCT using Bayesian inference method with the b γbld SR method in [7] and the additive NSCT method [5]. The
2 obtained images corresponding to each method are displayed xld − Eqk (y) [ybld ] +tr (cov k [ybld ])−1 , (27) in Fig. 4(d)-(f). p To assess the spatial improvement of the pansharpened images, we use the correlation of the high frequency compowhere ω = āω /c̄ω is the mean of the prior distribution nents (COR) which takes values between zero and one (the on the parameter ω and λω = āω /(āω + p/2), with ω ∈  higher the value the better the quality of the pansharpened ld ld ld αb , βb , γb . Eqs. (25)-(27) indicate that the inverse of image). Spectral fidelity was assessed by means of the peak the means of the parameters can be expressed as convex signal-to-noise ratio (PSNR), the Structural Similarity Incombinations of their prior values and their maximum likedex Measure (SSIM), an index ranging from −1 to +1, with lihood (ML) estimates and that λω , taking values in the +1 corresponding to exactly equal images, and the ERGAS interval [0, 1), can be understood as the confidence on the index, a global criterion for which the lower the value, the inverse of the mean of the prior distribution on the paramehigher the quality of the pansharpened image. See, for inters. So when λω are equal to zero, no confidence is placed stance, [1] for a detailed description of those measures. Taon the given values of the hyper-parameters and ML estible 1 shows the corresponding quantitative results. The promates are used, making the observations fully responsible of posed method with prior knowledge provides the best results the parameters estimation, while when they are asymptotifor each measure. The COR values reflect that all methods cally equal to one, the prior knowledge of the mean is fully are able to incorporate the details of the PAN image into enforced (i.e., no estimation of the hyper-parameters is perthe pansharpened one, although the SR method in [7], see formed). Since, as already said, the bands are related, we Figure 4(e), introduced less details in the band 3 (blue) since are going to study the relationships between the coefficients it contributes less to the PAN image and more into the band and use this relationships as prior values for the mean of 2 (green) since it has the highest contribution, which is rethe hyperparameters, with a high confidence, as will be deflected as a greenish color near the edges of the image. The tailed in the next section. Since the calculation of cov k [ybld ] NSCT method in [5] incorporates details in all the bands is very intense, we chose to approximate, only for parambut produces a noisy image, depicted in Figure 4(d). The eter estimation, W (uld )k in Eq. (24) by diagonal matrix, b p proposed method (Figure 4(f)) is able to incorporate details P p k k in all the bands while controlling the noise and avoiding the W (uld (uld b ) = 1/p b ) (i) × Ip . i=1 1/ 1133 (a) Original (b) Observed MS (c) Observed PAN (d) Add-NSCT in [5] (e) SR in [7] (f) Proposed method Figure 4: Results for the synthetic image (a) Observed MS (b) Observed PAN (c)Add-NSCT in [5] (d) SR in [7] (e) proposed method Figure 5: Results for the SPOT5 image color bleeding effect. The spectral fidelity measures shows that the proposed method performs better than the competing methods, which is also clear from the image in Figure 4(f). The PSNR for the proposed method is about 10dB higher than for the Add-NSCT method in [5] and from 2dB to almost 6dB higher than for the SR method in [7], with remarkable high SSIM and low ERGAS values which reflect the high quality of the resulting images. In a second experiment, the method was tested on a real SPOT5 dataset. Figure 5(a) shows a region of the RGB color image representing bands 1 to 3 of the MS image. Its corresponding PAN image is depicted in Figure 5(b). Visual inspection of the resulting images, displayed in Figures 5(c)(e), reveals similar conclusions to the obtained for the synthetic image. The proposed method calculates the parameters automatically and provides the best result, preserving the spectral properties of the MS image while incorporating the high frequencies from the panchromatic image and controlling the noise in the image. The NSCT method in [5] (Fig.5(c)) provides a detailed image but quite noisy, the SR method in [7] provides good details for bands 1 and 2, see Figure 5(e), but not for bands 3 and 4 since the PAN image does not cover those bands. This is why the blue color in Figure 5(e), seems to be vanished. The proposed method successfully preserves the colors, incorporates the details from the PAN image into the pansharpened image and controls the noise in the images. 6. CONCLUSION In this paper, a new pansharpening method based on superresolution reconstruction and non subsampled contourlet transform has been presented. The proposed method generalizes the fusion strategy of the panchromatic and multispectral images in contourlet based methods. The proposed fusion algorithm is based on the Bayesian modeling and incorporates a solid way to incorporate the details of 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 relationship between the contourlet coefficients has been examined. The use of this knowledge by the algorithm improved its performance. The method also estimate automatically the HRMS image and the unknown hyperparameters. The efficiency of pansharpening methods has been evaluated by means of visual and quantitative analysis, for synthetic and real data. The proposed method preserves the spectral properties of MS image while incorporating the high frequencies from the panchromatic image and controlling the noise in the image. Based on the presented experiments, the proposed method does significantly outperform NSCT-based and TV-based super-resolution methods. REFERENCES [1] I. Amro. Multispectral Image fusion using Multiscale and Super-resolution methods. PhD thesis, Dept. of Computer Science and Artificial Intelligence, Universidad de Granada, 2011. [2] A. L. da Cunha, J. Zhou, and M. N. Do. The nonsubsampled contourlet transform: Theory, design, and applications. IEEE Trans. Image Proc., 15(10):3089–3101, 2006. [3] I. Amro, J. Mateos, and M. Vega. General contourlet pansharpening method using Bayesian inference. In 2010 European Signal Processing Conference (EUSIPCO-2010), pages 294–298, 2010. [4] M. Song, X. Chen, and P. Guo. A fusion method for multispectral and panchromatic images based on HSI and contourlet transformation. In Proc. 10th Workshop on Image Analysis for Multimedia Interactive Services WIAMIS ’09, pages 77–80, 6–8 May 2009. [5] M. Lillo-Saavedra and C. Gonzalo. Multispectral images fusion by a joint multidirectional and multiresolution representation. Int. J. Remote Sens., 28:4065 – 4079, 2007. [6] D. D.-Y. Po and M. N. Do. Directional multiscale modeling of images using the contourlet transform. IEEE Trans. on Image Processing, 15(6):1610–1620, 2006. [7] 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. [8] N. Galatsanos. A majorization-minimization approach to total variation reconstruction of super-resolved images. In 16th European Signal Processing Conference (EUSIPCO 2008), August 2008. 1134