A Framelet-Based Iterative Pan-Sharpening Approach

Pan-sharpening is used to fuse multispectral images and panchromatic images to produce a multispectral image with high spatial resolution. In this paper, we design a new iterative method based on framelet for pan-sharpening. The proposed model takes advantage of the upsampled multispectral image and a linear relation between the panchromatic image and the latent high-resolution multispectral image. Since the sparsity of the pan-sharpened image under a B-spline framelet transform is assumed, we regularize the model by penalizing l1 norm of a framelet based term. The model is solved by a designed algorithm based on alternating direction method of multipliers (ADMM). For better performance, we propose an iterative strategy to pick up more spectral and spatial details. Experiments on four datasets demonstrate that the proposed method outperforms several existing pan-sharpening methods.


Introduction
An optical satellite usually acquires two images describing the same scene almost simultaneously, which are called multispectral (MS) image and panchromatic (PAN) image respectively.The former is a multichannel image with low spatial resolution, while the latter is a single channel image with rich spatial details.Despite the fact that MS image can be of over eight bands and the resolution of PAN image can be less than half a meter, their superiority cannot be synthesized in one image due to physical and technological constraints.However, pan-sharpening techniques are capable of creating a multichannel image with high spatial resolution out of these two images, which is of great importance for remote sensing.More specifically, pan-sharpening plays an important role in the interpretation of remote sensing scenes and can be used as a preliminary step of various remote sensing tasks such as object recognition [1], change detection [2] and so on.Therefore, these techniques attract much attention of scientific community.
Among various pan-sharpening methods having been proposed in the literature, many of them can be put into two main categories: component substitution (CS) and multiresolution analysis (MRA).CS mainly involves three steps: firstly a spectral transformation of MS image, then a replacement of its spatial component with PAN image, and finally an inverse transformation.This class includes classical methods such as intensity-hue-saturation (IHS) [3], principal component analysis (PCA) [4], Gram-Schmidt spectral sharpening (GS) [5] and recent methods based on mean information [6] and image matting model [7].As for MRA, it focuses on an injection process in which spatial details extracted from PAN image are added to upsampled MS image.Examples of this class include wavelet transform based methods [8][9][10], Laplacian Pyramid based methods [11], and methods based on other transforms [12,13].
Different from the two categories mentioned above, there are also other types of pan-sharpening methods.This family includes those based on Bayesian paradigm [14], total variation [15,16], gradient operator [17,18], sparse representation [19], super-resolution techniques [20], convolution neural network [21] and so on.Recently, Deng et al. [22] propose a novel variational model for pan-sharpening in which intensity function of the unknown image is considered from a continuous point of view.The related continuous function is made up of two components.Assumption is made that the former lies in a Reproducing Kernel Hilbert Space (RKHS) while the latter can be approximated by linear combination of approximated Heaviside functions (AHF).This model outperforms several state-of-the-art pan-sharpening methods according to experiments on two datasets.However, its good performance relies on a large amount of computation, which results in rather long running time.
In this paper, a new iterative algorithm for pan-sharpening is proposed as an attempt to simplify RKHS method.We make use of the information from MS image by generating its upsampled form.The linear relation [22,23] between PAN image and bands of the image to be estimated is also considered in the model.A framelet based term is introduced as regularization.Besides, we adopt an iterative strategy similar to [22] to improve performance of the algorithm.The framework of the proposed approach can be seen in Figure 1.By utilizing real data from Pléiades, Quickbird, WorldView-2 and SPOT-6, we compare several existing pan-sharpening methods with the proposed method.These results show spatial and spectral fidelity of the proposed method.Meanwhile, they confirm a much less running time than RKHS method, suggesting that framelet based regularization is more mature and easier to compute.The rest of this paper is arranged as follows.Section 2 reviews the related work [22].Then in Section 3 the new iterative algorithm is presented.Section 4 is a display of visual and numerical experimental results together with some discussions about the proposed method.Finally, conclusion is drawn in Section 5.

Related Work
To begin with, we introduce several notations to be used throughout this paper.Let MS ∈ R m 1 ×n 1 ×T be the original multispectral image with T bands, each band denoted as MS i ∈ R m 1 ×n 1 .Let MS ∈ R m 2 ×n 2 ×T be the upsampled multispectral image, each band denoted as MS i ∈ R m 2 ×n 2 .MS represents the high-resolution multispectral image to be estimated, with MS i ∈ R m 2 ×n 2 being the ith band of it.Moreover, the original panchromatic image is P ∈ R m 2 ×n 2 .
Deng et al. propose a new iterative pan-sharpening algorithm [22], which is an extension of their previous work on super-resolution [24].They view the pan-sharpening problem as an intensity estimation process for the unknown image MS.The intensity of MS is modelled as a hidden continuous function.It consists of two different components, a smooth one and a non-smooth one.The former is assumed to be an element of a special function space called Reproducible Kernel Hilbert Space (RKHS), while the latter is formulated as a linear combination of approximated Heaviside functions (AHF).
Specifically, let f i , i = 1, 2, ..., T, be the underlying continuous intensity function corresponding to the ith band of MS.Without loss of generation, the domain is restricted such that . These two series lie in two different RKHS, with φ ν,i (z), ν = 1, 2, ..., M, and ξ s,i (z), s = 1, 2, ..., n, being basis functions in each RKHS respectively.d ν,i , ν = 1, 2, ..., M, together with c s,i , s = 1, 2, ..., n, are the corresponding coefficients.As for the non-smooth component, a family of Heaviside functions is considered to implement the approximation.However, since Heaviside functions are singular at the origin, which makes differentiation impossible to be taken there, they are approximated by a smooth form, i.e., approximated Heaviside functions (AHF).In 2D case, this smooth alternative takes the form of ψ((cosθ j,i , sinθ j,i ) • z + c ρ,i ), actually representing an edge with θ elevation at a location specified by c ρ,i .It is a generalization of the 1D form ψ(x) = 1 2 + 1 π arctan( x ξ ), where ξ is a positive parameter to control smoothness.Therefore, f i , i = 1, 2, ..., T, is finally modelled as: For more details, see [22,24].By evaluating on a fine grid, (1) can be discretized so that it becomes a form involving simple matrix multiplication.
matrices whose entries are generated by evaluating φ ν,i (z), ξ s,i (z), and ψ((cosθ j,i , sinθ j,i ) • z + c ρ,i ) on a fine grid.Then each band of the desired multichannel image with high spatial resolution can be computed as follows: Hence the pan-sharpening problem is converted into a problem of coefficient estimation.In [22], the coefficients are computed by minimizing the following model: where N is the number of pixels that MS contains.µ, λ 1 , and λ 2 are positive parameters.MS i is generated by an upsampling process via GS method [5].
, the discretization process is done on a coarser grid.ω i are weights reflecting the contribution of each band of MS to the linear combination which approximates the panchromatic image P.This linear approximation and its variants are assumed by many methods, e.g., [23,25].
After the coefficients d i , c i , β i are obtained, the high-resolution multispectral image can be computed by (2).Furthermore, model (3) is combined with an iterative strategy.It will be detailed in Section 3 since it is used in the method proposed in this paper as well.

The Proposed Method
Empirical results in [22] show that RKHS method outperforms several state-of-the-art pan-sharpening methods both in the perspective of spatial and spectral fidelity.However, the RKHS based model of MS is quite complicated, since it is not easy to implement in actual calculation.In addition, the regularization terms of (3) involve not only l 1 norm of β i but also quadratic function of c i (both summed by band), which adds to computational burden.As a result, the algorithm is rather time consuming.Therefore, it is natural to raise a question.Can we simplify model (3) without loss of its advantage, i.e., high-level spatial and spectral performance?A possible way is not to make efforts to build a complicated model of MS, but seek for another technique of regularization instead.This is exactly why we turn to framelet.
Piecewise smooth functions, for instance images, can be sparsely approximated by framelet system efficiently [26].As a result, framelet techniques are used in literature to address the problems like image restoration (e.g., [26][27][28]).Recently, it is also applied to pan-sharpening.For instance, a framelet based MRA scheme is considered in [29].A variational model [30] based on assumptions related to framelet coefficients, geometry keeping, spectral preserving and the sparsity of the image in the framelet domain is also proposed.They are able to obtain good results.Nevertheless, in this section, we consider a combination of variational model and iterative strategy instead of building a complicated model, i.e., we try to build a simple framelet based variational model and then combine it with an iterative strategy to yield a novel and effective approach for pan-sharpening.
In discrete case, let W and W T denote fast framelet decomposition and construction respectively.They are constructed by unitary extension principle (UEP) [31], and satisfy the relation W T W = I.In this paper, we use the piecewise linear B-spline framelets [26].An L-level framelet transform of an image u [26] can be denoted as: where W l,j u denotes the coefficients of u in framelet band j at level l under the framelet transform, and I is the index set of all framelet bands.Throughout this paper, we empirically set L = 1.For more theoretical details of framelet, see e.g., [32].Consequently, by employing framelet based regularization, (3) can be modified.Each band of MS is computed by minimizing the following function: where α and λ i are parameters.Dot product is used in the third term since there are more than one framelet band, as (4) suggested.Coefficients for bands of MS, i.e., ω i , are estimated automatically by a linear regression [23] between the original multispectral image and downsampled panchromatic image.The upsampled multispectral image MS i is generated via GS method.Note that this model can also be viewed as an extension of the so-called analysis based model for image restoration (e.g., [33,34]).Each l 1 term of ( 5) is in accordance with that defined in the analysis based model.Concretely, it can be expressed as This expression is the case where 1-level framelet transform is imposed on MS i , as we emphasized above.
Model (5) can be solved by methods such as primal-dual method [35] and ADMM [36] efficiently.We choose ADMM here, whose application covers a wide range of image processing, such as image denoising [37], image super-resolution [38], tensor completion [39], image destriping [40,41] and so on.Its convergence is guaranteed by many works such as [42,43].Due to non-smoothness caused by l 1 term, the first step is to rewrite (5) as an equivalent form through substitution of variables: Then we can obtain the augmented Lagrangian of (7), i.e.,

L( MS
where D i and E i are Lagrangian multipliers, β 1 and β 2 are two positive parameters.Now we denote According to ADMM, problem (7) can be solved by implementing the following iterative scheme: (1) For i = 1, ..., T, update each u (k) i by solving: (2) For i = 1, ..., T, update each V (k) i by solving: ( by solving: ) ).
Note that ( 9)-( 11) have closed-form solutions.Using soft-thresholding operator (e.g., [44]) T τ , u (k+1) i can be rewritten as: where T τ (ν) is defined entry-wise by As for (10) and (11), they can be solved easily as follows: In order to facilitate illustration of the proposed algorithm, we summarize these steps of ADMM as Algorithm 1.Note that for simplicity of form, we write each iteration of ADMM in Algorithm 1 in an order slightly different from what we mentioned above, but it is easy to validate that they give the same results.

end for end while
Although Algorithm 1 itself is a complete algorithm, there is still room for improvement.Thus an iterative stategy is considered.For accordance of notations, let P (1) = P and MS (1) = MS.Similarly, denote the first output of Algorithm 1 as I (1) .After obtaining MS, we compute In addition, let MS (2) = U(MS-D( MS)).D represents a downsampling operator, and U represents the upsampling process by which MS is generated.Now we view P (2) and MS as new inputs for Algorithm 1 instead of the original P and MS.The resulting new output can be denoted as I (2) .By repeating this strategy, we obtain a series of I (j) , j = 2, 3, ..., γ.Then the sum of all I (j) , j = 1, 2, ..., γ, is taken as the final high-resolution multichannel image.This iterative strategy is adopted in not only pan-sharpening [22] but also image super-resolution [24].Now we can summarize the procedures above as Algorithm 2: Algorithm 2 The proposed iterative pan-sharpening algorithm Input: panchromatic image P, multispectral image MS, ω i α, λ i , β 1 , β 2 .
Output: high-resolution multispectral image MS.
(3) Update P (j) by P (j+1) = P Note that γ is the number of outer iterations.The downsampling process D in each iteration is completed through a combination of two steps: compute the modulation transfer function of I (j) with Gaussian filter and then interpolate it to the size of MS in a "nearest" way [22].And the upsampling process is done by GS method.

Results and Discussion
In this section, we firstly utilize four datasets to compare the proposed method with several pan-sharpening methods.After that, discussions related to the number of outer iterations of Algorithm 2 are discussed.Results on time cost are presented as well.
Since high resolution multispectral images are not available in the datasets, we follow Wald's protocol [46].Therefore, the original multispectral images in the datasets are treated as ground truth.The scale ratio is 4, thus the simulated low-resolution multispectral images (4 bands) are of the size 128 × 128, 256 × 256, 200 × 200, and 256 × 256 respectively.Each of them is downsampled from the corresponding ground truth in the same way as that in Algorithm 2, i.e., filter the ground truth by a Gaussian filter matched with the modulation transfer function (MTF) and then downscale it by "nearest" interpolation.As for each P, it is generated by combining bands of the ground truth linearly.
Parameters of the proposed algorithm are empirically set as follows.We use 1-level piecewise linear B-spline framelet.For each band, λ i is equally set as: Other model parameters in (8) are set as α = 1.5, β 1 = 0.5, and β 2 = 0.5, with coefficients ω i estimated by linear regression.In addition, the number of the external iterations in Algorithm 2, i.e., γ, is set as 5.For different datasets, these settings may not always be the best choice, but we unify them to display stability of the proposed method and also to save efforts of tuning.
Methods compared with the proposed method comprise some classical pan-sharpening methods (PCA [4], GS [5], high-pass filtering (HPF) [47], and modulation transfer function-generalized Laplacian Pyramid (MTFGLP) [48,49]) and different kinds of recent state-of-the-art methods (pan-sharpening with hyper-Laplacian prior (PHLP) [50], nonlinear intensity-hue-saturation (NIHS) [51], and RKHS [22]).All the experiments are conducted in MATLAB on a laptap with 4GB RAM and 1.70 GHz Intel(R) Core(TM) i5-4210U CPU.It is obvious that all of the other methods perform better than PCA method both in terms of spatial quality and spectral fidelity according to Figures 2-4.GS method outperforms PCA method significantly but generally fails to avoid great spectral distortion, which is most visible on Quickbird dataset and Pléiades dataset.NIHS method and PHLP method preserve spectral characteristics quite well.However, from the perspective of spatial details, they tend to generate excessive smooth results, thus the pan-sharpened images provided by them lack much sharp spatial information compared with the reference images.

Visual Comparison
HPF method and MTFGLP method are visually without much spectral distortion and are able to keep more spatial details.However, a closer look at their resulting images, e.g., in Figure 3, suggests that they are not able to provide as many spatial details as the proposed method does.Finally, it is noticeable that RKHS method and the proposed method achieve the best visual performance on the last three experimented datasets.Actually, visual comparisons show little disparity for these four methods.However, we will demonstrate that the proposed method gives the best quantitative results in Section 4.2.

Quantitative Comparison
Several quantitative indices are employed to report the performance of different pan-sharpening methods.To evaluate spectral distortion, we use spectral angle mapper (SAM) [52], erreur relative globale adimensionnelle de synthése (ERGAS) [52], universal image quality index (Q) [53] together with its vector extension Q4 [54], and relative average spectral error (RASE) [55] (the larger Q and Q4 and the smaller SAM and ERGAS, the better performance).Correlation coefficient (CC) [56], with 1 being its ideal value, acts as spatial quality metric.Meanwhile, we use peak signal-to-noise ratio (PSNR) and root mean square error (RMSE) as metrics of fusion accuracy.Generally speaking, better performance is achieved when PSNR is larger and RMSE is smaller.
In each experiment, most of the compared methods require an upsampled multispectral image as an input.Unless specially specified in the literature, we unifiedly generate them by interpolating via a kernel function which is a polynomial with 23 coefficients [57].
The quantitative results of four datasets with regard to eight metrics are reported in Tables 1-4.They clearly validate that the rest of the methods outperform PCA method as visual comparison in Section 4.1 preliminarily confirmed.GS method maintain spectral fidelity well on Quickbird dataset and WorldView-2 dataset, while on the rest two datasets it performs not so well.Compared with methods such as MTFGLP method and HPF method, PHLP method and NIHS method give comparable results with respect to SAM.However, there is still a gap when it comes to metrics reflecting spatial quality and fusion accuracy.Similar to GS method, HPF method is unable to preserve enough characteristics on Pléiades dasaset and SPOT-6 dataset from a spectral point of view.If not taking the proposed method into consideration, RKHS method or MTFGLP method performs the best.However, the proposed method consistently achieves better quantitative performance than them in terms of all metrics (except for SAM) on the four tested datasets.These observations are sufficient to demonstrate that the proposed method preserves spectral information and sharp spatial details accurately.In Algorithm 2, we use an iterative strategy to improve the performance of Algorithm 1.The number of this outer iterations, i.e., γ, is set to 5. It is meaningful to inspect how the performance of Algorithm 2 changes as γ increases.
In Figure 6, we present I (j) , j = 1, 2, 3, 4, 5, each being the result of the jth outer iteration of Algorithm 2 with respect to Pléiades dataset.To focus on spatial details, only the first channel of each image is shown.For each of the last four images, we subtract its intensities by the corresponding smallest value in the channel before normalization, since these images contain negative intensities which cannot be plotted by MATLAB directly.For better visualization, the normalization of these four images are implemented by dividing the largest intensity of I (2) since their absolute values are rather small compared with I (1) .From these visual results, we know that the outer iteration in Algorithm 2 is able to pick up more image details.
From Tables 5-8, quantitative results of the proposed method with γ changing from 1 to 6 are listed.When γ = 1, the proposed method essentially becomes Algorithm 1, which makes it convenient to compare it with the performance of Algorithm 2 directly.As expected, after combining Algorithm 1 with the iterative strategy, better quantitative performance is achieved, which can also be inferred from Figure 6.    3) ; (e) the first channel of I (4) ; (f) the first channel of I (5) .
We observe from Tables 5-8 that the best quantitative performance is achieved at γ = 5 when γ is not too large, i.e., smaller than 6.An exception can be noticed in the case of SPOT-6 dataset, where the best performance is observed at γ = 4.However, for the rest of the datasets, the case of γ = 5 outperforms other cases.
When we let γ goes even larger, things become different.This can be seen in Figure 7, which shows how indices vary when γ goes from 1 to 20.Since situation for all the tested datasets and all the indices is similar, we only present two indices related to Pléiades dataset here.From Figure 7, it is obvious that the curves of indices fluctuate when γ is larger than 3.However, it is also noticeable that there is not much improvement for the performance of Algorithm 2 corresponds to larger γ compared with the case of γ = 5.To explain this phenomenon, an iterative regularization algorithm [58] may helpful, whose output signals become noisy if the outer iterations are not stopped properly, i.e., smaller than a threshold.And the authors explain that the phenomenon of noisiness is due to the influence of noisy input image.Despite the fact that the output images of Algorithm 2 do not necessarily degrade all the way when γ becomes larger than 5, the explanation of the similar phenomenon in [58] is enlightening.Therefore, we suggest that the output of Algorithm 2 is influenced significantly by the blurred upsampling multispectral image when γ is large enough, which leads to the fluctuation of numerical indices.Since the success of our iterative strategy has been preliminarily demonstrated by numerical results, an analysis of this conjecture can be explored in our future research.We simply set γ as 5 as a tradeoff of computational burden and the performance of the proposed algorithm, and also as a strategy to save the efforts of tuning.

Time Comparison with RKHS Method
Since we mentioned in Section 3 that the proposed model is supposed to simplify model (3), it is reasonable to expect a reduction in running time when comparing these two methods.Time cost of both methods for the four datasets can be examined in Table 9, whose values are presented in the unit of second.We point out that these results are obtained by running on the same laptap with 4GB RAM and 1.70GHz Intel(R) Core(TM) i5-4210U CPU as mentioned in the beginning of Section 4. From the table, it can be validated that the proposed method is much less time consuming than RKHS method, which means the goal of reduction in computational burden and running time is achieved.

Conclusions
In this paper, a new iterative pan-sharpening approach is proposed for the fusion of panchromatic image and multispectral image.The proposed variational model inherits the framework of RKHS method but is essentially different.Instead of paying attention to modelling the unknown high resolution multispectral image in a complicated way, we utilize framelet technique in image restoration for more effectiveness of regularization.An iterative scheme similar to [22] is also employed to improve performance of the proposed algorithm.Experiments on data from Quickbird, Pléiades, WorldView-2, and SPOT are conducted for visual and numerical assessment.The results demonstrate that the proposed method outperforms several state-of-the-art pan-sharpening methods both visually and quantitatively.Meanwhile, it succeeds in reducing time cost compared with RKHS method.
image and panchromatic image respectively.• is defined by (5) in the paper.• U and D are upsampling and downsampling process respectively.• is the number of external iterations.

Figure 1 .
Figure 1.Framework of the proposed iterative pan-sharpening approach.

Figures 2 -
Figures 2-5 show visual results obtained by conducting experiments on four datasets aforementioned.Each set of figures contains output images produced by eight different methods.Each ground truth is also presented as reference.For better visualization, we show local enlarged images at the bottom-left corner of each output image.It is obvious that all of the other methods perform better than PCA method both in terms of spatial quality and spectral fidelity according to Figures2-4.GS method outperforms PCA method significantly but generally fails to avoid great spectral distortion, which is most visible on Quickbird dataset and Pléiades dataset.NIHS method and PHLP method preserve spectral characteristics quite well.However, from the perspective of spatial details, they tend to generate excessive smooth results, thus the pan-sharpened images provided by them lack much sharp spatial information compared with the reference images.HPF method and MTFGLP method are visually without much spectral distortion and are able to keep more spatial details.However, a closer look at their resulting images, e.g., in Figure3, suggests that they are not able to provide as many spatial details as the proposed method does.Finally, it is noticeable that RKHS method and the proposed method achieve the best visual performance on the last three experimented datasets.Actually, visual comparisons show little disparity for these four methods.However, we will demonstrate that the proposed method gives the best quantitative results in Section 4.2.

Figure 7 .
Figure 7. Performance of Algorithm 2 as γ increases, represented by SAM (red) and ERGAS (blue) with respect to Pléiades dataset.

Table 1 .
Quantitative results for Quickbird dataset.

Table 5 .
Quantitative results for Quickbird dataset with different numbers of outer iterations.

Table 6 .
Quantitative results for Pléiades dataset with different numbers of outer iterations.

Table 7 .
Quantitative results for WorldView-2 dataset with different numbers of outer iterations.

Table 8 .
Quantitative results for SPOT-6 dataset with different numbers of outer iterations.

Table 9 .
Running time(s) of RKHS and the proposed method on four datasets.