An Innovative Pansharpening Method Based on MRF Strategy

An innovative pansharpening method is proposed to selectively extract more useful information from the original images to produce a new image with higher resolution information. The standard PCA is employed as a decorrelation tool to separate the spectral and spatial information in MS images. In order to reduce the spectral distortion of fused image, we decompose the first principal component (PC1) of multispectral (MS) images and panchromatic (PAN) images using nonsubsample shearlet transform (NSST) to achieve effective detailed information; a novel energy function, including the interand intrainformation between subbands, has been established to take full account of the local dissimilarity between MS and PAN images, and the reasonable coefficients are selectively chosen based onMarkov random field (MRF). It is found that the simulated image by the newmethod is more close to the real image and more clear and with more detailed information compared with other popular methods reported recently, which means that our new method can effectively improve the efficiency and quality during the fusion image process.


Introduction
Image fusion is an important part of information fusion that synthesizes two or more source images obtained from different sensors to a fused one of greater quality [1].A large collection of pansharpening methods have been designed and developed over the last two decades [2].These methods are mainly divided into three categories: projection-substitution (the intensity hue saturation (IHS) [3] and the principal component analysis (PCA)), relative spectral contribution, and the multiscale resolution analysis concept [4].The Gram-Schmidt (GS) method [5], the principal component analysis (PCA) [6] method, and the generalized intensity hue saturation (GIHS) [7] method serve as some of the most popular projection-substitution based methods.Besides, Brovey [8] and the additive wavelet luminance proportional (AWLP) [9] methods also contribute to the relative spectral and ARSIS analyses, respectively.In addition, recently, two novel pansharpening methods [10,11] which are elaborately designed based on sparse representation have been proposed with quite satisfactory performance.
Being simple in calculation and good in performance, the projection-substitution based methods have been adopted in many cases.The PCA method is their representative, and in the method the first principal component (PC1) with the largest variance is replaced by the histogram-matched PAN image.However, the high resolution PAN image is not ideally equal to the substituted component from the corresponding low resolution MS images, which may usually lead to spectral distortion.To improve the performance of PCA and make more use of the consistency of the images, the theory of Markov random field (MRF) is used to explore the correlations of the local pixels.The MRF is an ideal tool that is able to incorporate much information provided by the neighborhood because of its flexible energy function.Therefore, by establishing a novel energy function which includes the inter-and intrainformation between pixels, we propose a generalized fusion method based on MRF, which extends the MRFFCM [12] to domain of pansharpening.Besides, in order to get the directional and detailed information of both MS and PAN images, NSST is used as a multiresolution analysis tool to decompose the corresponding images.Because of the property of shift-invariance, NSST can effectively reduce the Gibbs phenomenon and improve computational accuracy [13,14].
This paper focuses on a particular application task of the synthesis of remote sensing images coming from the optical satellites such as IKONOS, QuickBird, and WorldView-2.This is called the pansharpening, that is, combining the MS images which have characteristics of low spatial resolution but high spectral resolution and PAN image which has reverse properties.The purpose of pansharpening is to obtain a fused image as close as possible to the image captured by the multispectral sensor if it would work at the high spatial resolution of the PAN.
In this paper, we propose a novel pansharpening method in which a robust selection strategy based on MRF is included.In the method, MRFFCM is employed between high frequency information of PC1 originating from the low resolution MS images and the high frequency of PAN acquired after NSST decomposition.Thus, a new fused PC1 (PC1 fused ) is constructed.MRFFCM serves as a tool to generate two membership matrices of PC1 and PAN which the selection strategy is based on, and the novel energy function is applied to MRFFCM to obtain the weighting planes by discovering the relationship of the two matrices.Due to the fact that the local spatial information between PAN and MS images is involved, PC1 fused can be constructed with much deliberateness.Hence, the proposed method can improve quality of the fused MS images both qualitatively and quantitatively.This procedure takes the local dissimilarities between PAN and MS images into consideration and significantly improves the phenomena of spectral distortion.
The rest of the paper is organized as follows.Section 2 reviews the theory of the NSST and MRF briefly and then describes the novel energy function in the MRF and the main procedure of the proposed method.The experimental results, the real date set, and the performance of the corresponding method are discussed in Section 3. Finally, conclusions are drawn in Section 4.

Methodology
The MRF-based pansharpening method comprises several parts: (i) standard PCA and nonsubsample shearlet transform (NSST), (ii) analysis of high frequency of first principal component (PC1  ) and high frequency subband of PAN obtained from part (i) based on MRF strategy, and (iii) proposing the new energy function and applying MRFFCM to the corresponding bands.A flow chat summarizing the whole procedure is provided in Figure 3.The preparation work is that the original MS images and PAN image need first to be registered with each other and then the registered MS images are resampled to the same scale of PAN images.

NSST.
The discrete shearlet transform can be obtained by subsampling the shearlet on a proper discrete set.NSST consists of two parts, and the image is first decomposed by the nonsubsampled pyramid (NSP), which produces one lowfrequency subband and +1 high frequency subbands whose sizes are all the same as source image;  denotes the number  of decomposition levels.NSP ensures the low-frequency component available iteratively to capture the singularities in the image.After that, the high frequency subbands are decomposed by the shearing filters with  stage at each scale, which will yield 2  directional subbands with the same size as the source images.In frequency domain, each shearlet is supported on a pair of trapezoids, of approximate size 2 2 × 2  , oriented along lines of slope 2 − .Figure 1 shows the proceeding of NSST, and the frequency partition of NSST is shown in Figure 2.

MRF.
Traditional fusion strategies are fused so locally that regions of structures may not be consistent.The MRF is a useful tool for analyzing and integrating contextual information.The fusion strategies of injecting the high resolution information are explored by the modified MRFFCM model.The method is more suitable for the purpose of human visual perception and computer processing due to its consideration of local information rather than arbitrary pixels.
Let  = { = (, ) | 1 ≤  ≤ , 1 ≤  ≤ , , , ,  ∈  + } be the set of nodes corresponding to the pixel site, where  and  denote the height and width of the image, respectively.According to the MRF theory, relationship between sites in  are constrained in a neighborhood system denoted by  = {  ,  ∈ }, where   is a set of pixel sites in  neighboring .Let ( = ) denote the a priori probability of label ; according to the Hammersley-Clifford theorem, a Gibbs distribution of MRF  is given by where  is a normalizing constant (the partition function) given by where () is the energy function comprised of some clique potentials.In order to extend the MRFFCM to the domain of pansharpening, it is important to establish a proper energy function in a certain scenario.The novel energy function will be described in the next section.

The Novel Energy Function.
The MRF-based fusion model enables the accurate estimation of fusion weights for the detailed information of PC1 and PAN.In order to make full use of the relationship of the spatial domain, a new energy function under the concept of Markovian neighborhood is described as follows: where  intra and  inter denote intra-and inter-subband weights and PAN   represents the th high frequency subband of PAN image.DI is short for the difference of images and is given as As DI denotes the high value of the corresponding pixels, these pixels always represent the edges and detailed information of the images.Thus the new energy function cannot only describe the intrarelationship of the high frequency bands, but also contain the interinformation through the DI procedure, and therefore the modified MRF can explore the links of detailed information between PC1   and PAN   .Then after FCM procedure to each high frequency subband of PC1   and PAN   , we can get weighting planes as  1 ,  2 , . . .,   .These planes will decide the detailed information of the fused MS.
In order to utilize the information of the neighborhood system adequately, we view the mean-like field as an approximation of the MRF as where   = mean ∈  {  }.

MRFFCM Based Fusion Strategies.
We first apply the standard PCA transform to MS images to acquire PC1 and then decompose PC1 and PAN image by NSST to obtain scale-by-scale subbands as follows: where PAN  and PC1  denote the low spatial frequency subband of PAN and PC1, respectively.In the same way, PAN   and PC1   denote the th high spatial frequency subband of PAN and PC1.
Finally, FCM is used to obtain the weighting planes of high frequency subband of PAN image.In fact, this procedure is applied to achieve the maximum a posteriori (MAP) estimation for obtaining the weighting planes   and chooses the coefficients which are more intense and more similar to PC1  .MAP can be described as where () is calculated by MRF in (1); let  denote the th iteration and the conditional probability ( | ) is described by Gaussian mixture model (GMM) as follows: where    and    are the membership degree and the standard deviation of the corresponding classes.   denotes the distance matrix; in addition,  = ,  represents different classes, and according to (2), (7), and (8), we can formulate the new membership degree as follows: And then we can utilize FCM to calculate the objective function    as follows: Update the mean membership  +1  and standard deviation  +1  , respectively, as given in the following: The iteration will stop until the objective function    converges to some extent and output {  }.Let  denote the threshold of the relationship matrix {  }; the selection strategy of weighting planes   is defined as where   =   means the class which possesses the more intense energy and presents more similar to PC1  is chosen.And then the high frequency PAN   can be selectively injected to the PC1 with the weighting planes as follows: According to ( 14), PC1 fused is used instead of PC1 and the inverse PCA transform is performed to obtain the new fused image.
The procedure of the proposed method can be summarized as follows.
Step 1. Coregister the original low resolution MS images with the high resolution PAN image and resample the MS images to the same resolution of PAN image.
Step 2. Apply the standard PCA to the resampled MS images to acquire the PC1 image.
Step 3. Apply NSST to PC1 and PAN image to obtain PC1  and the detailed information PAN   ( = 1, 2, . . ., ), where  denotes the th subband of PAN.
Step 4. Employ MRFFCM on both PC1  and PAN   to construct PAN  fused .
Step 6. Make inverse PCA to acquire the fused MS images.

Experimental Results and Discussion
Based on three datasets originating from QuickBird and IKONOS satellites, this section evaluates performance of the proposed method.The two QuickBird datasets with 2.4 m resolution in MS and 0.6 m resolution in PAN are two scenes of Sundarbans, India, and they are referred to as Sundarbans1 dataset and Sundarbans2 dataset, respectively.The IKONOS data show a scene of Sichuan, China, with the spatial resolution for IKONOS data 4.0 m in MS and 1.0 m in PAN, and it is named as the Sichuan dataset.The sizes for all the three datasets are 256 × 256 × 4 in MS and 1024 × 1024 in PAN.
Adopted by Walds protocol [15], in the experiment, we first decimate the MS and PAN data by a factor of 4 using bicubic interpolation.Then, these decimated data are sharpened to acquire the fused 256 × 256 × 4 MS images.Finally, the fused MS images are evaluated with the original MS images which are treated as the reference high resolution MS images.We also use GIHS, Brovey, GS1, PCA, and AWLP as the comparison algorithms.The spectral angle mapper (SAM), the relative dimensionless global error in synthesis (ERGAS), and the quaternion theory-based quality index (Q4) [16] are used to quantitatively evaluate the performance of the algorithms.For the three quantitative indexes, SAM is used to measure the global spectral distortion by degrees and ERGAS to measure the overall spectral quality of fused MS images.Q4 measures both the spatial and spectral qualities of the fused MS images.SAM and ERGAS reflect the high quality with lower value, while the high value of Q4 represents good qualities of both spatial and spectral qualities of the fused MS images.
Figures 4-6 show the three datasets and their corresponding experimental results.roads and buildings.This is because after PCA is used to separate the spectral information, we improve the resolution of the resampled MS image by injecting high frequency information of PAN.In addition, MRFFCM is adopted to exploit the local spatial information between PC1  and PAN   .When observing the suburban areas results shown in Figure 5, we can see that the methods (c)-(h) combine Figures 5(a Finally, when it comes to the rural areas shown in Figure 6, the textures in the image also exist in a clear manner, and the similar conclusion can be drawn that spectral characteristics in Figure 6(h) are also compatible with Figure 6(b).
For the quantitative indexes of Figures 4-6, they are given in Tables 1-3, respectively.In the three tables, the last column "Ideal" denotes the ideal value for the three quantitative indexes.For each index, the best one among all the methods is highlighted in bold.In general, the proposed method leads to much better performance than others, and this is due to the fact that the novel energy function in MRFFCM involves an elaborate local relationship.It helps generate a more accurate membership and thus establishes a robust strategy.From both qualitative and quantitative results, it is demonstrated that the proposed method can generate a good fused MS image with high performance of clearer object features and low spectral distortion.

Conclusion
This paper introduced an innovative pansharpening method of applying MRFFCM to consider the condition of local dissimilarities.By establishing the novel energy function, the modified MRFFCM enables the accurate estimation of fusion weights for the salient and intense coefficients.
Compared with the conventional method, our new method could preserve the contour structures of the original images and deserves less spectral distortion.
In the future work, it would also be possible to explore the relationship between coarse subband and fine subband in terms of Markov theory and establish the link to guide the algorithms with more accuracy.
) and 5(b) in some extent.But it can be concluded that some spectral distortions present in Figures 5(c), 5(d), and 5(g) compared with those in Figure 5(b).The roads of forests in Figure 5(h) are more suitable for human vision system, as the roads and forest of Figure 5(h) are more distinct and smooth than those of Figures 5(e) and 5(f).

Table 1 :
Quantitative indexes for fusion results from Sundarbans1 dataset.

Table 2 :
Quantitative indexes for fusion results from Sundarbans2 dataset.

Table 3 :
Quantitative indexes for fusion results from Sichuan dataset.