Application of Morphological Component Analysis to Optical Image Fusion

: The image fusion technique is widely used in remote sensing. Its purpose is to provide comprehensive information without arte facts by combining the partial information from different source images. In this study, we propose a new model of images fusion with very high spatial resolution. We use the separation capacities of the Morphological Component Analysis (MCA) to extract the smooth and texture components of our images. These morphological components are then fused separately using the decomposition in the Laplacian pyramids for the smooth part and bivariate Hahn polynomials for texture part. Finally the image fusion is obtained through linear combination of merged smooth and texture components. The experiments carried out on IKONOS, LANDSAT and Quick Bird remote sensing images show the good performances of our method which has been compared to conventional methods. The performances obtained in our experiments are characterized by a small global metric such as ERGAS equals to 3.88 for IKONOS image and 3.65 for QuickBird image compared to 8.70 for IKONOS image and 6.97 for QuickBird for conventional HIS algorithms. We also have a mean loss of 15% for spectral information com pare to those of the conventional methods which revolve around 25%. The degradation of spatial information in order of 17% in contrast to conventional HIS algorithms which oscillate around 21%.


Introduction
The data fusion techniques permit the extraction and the combination of the information from multiple heterogeneous sources to assess the targeted phenomenon to which there is interest. They help greatly to make a decision. The application fields of the data fusion are varied and diverse: Medical imaging (Magtibay et al., 2016;Adali et al., 2015), economy (Barenboim and Pearl, 2016), information theory (Gagolewski, 2016), image processing (Paris and Bruzzone, 2015), etc. In remote sensing where the nature and the resolution of sensors are various and different, methods of image fusion implement several types of images: Panchromatic (PAN), Multispectral (MS), Hyperspectral (HS) and Radar (SAR). These methods combine either panchromatic and SAR images (Joshi et al., 2016), or panchromatic and HS images (Loncan et al., 2015) or panchromatic and multi spectral images (Vivone et al., 2014).
The fusion of the panchromatic and MS image known as pansharpening has been the subject of numerous publications in the literature (Joshi et al., 2016;Pandit and Bhiwani, 2015). The proposed methods can be classified into methods with injection of high spatial frequencies, methods based on substitution of components and methods based on a multiresolution analysis.
The first methods extract the details of PAN image by a high frequency filter and inject them into the MS image through a subtraction operation (Ranchin et al., 2003). The second methods submit the MS image to a transformation and substitute the panchromatic image to one of the components of the MS images before forming the fused image through the inverse transformation. The most common transformations are IHS (Intensity-Hue-Saturation) transformation for RGB color image (Rahmani et al., 2010;Lari and Yazdi, 2016), Principal Components Analysis (PCA) for multispectral images (Dou et al. 2007) and Gram-Schmidt method (Aiazzi et al., 2007) which is a generalization of the PCA. Finally, the third group of methods is based on a Multi Resolution Analysis (MRA) by Laplacian pyramids (Aiazzi et al., 2012) or by Wavelet, curvelets and contourlets Transforms (Lari and Yazdi, 2016;Otazu et al., 2005;Amolins et al. 2007). In these latter methods, multiresolution analysis is conducted on the panchromatic and MS images up to a given level. The different types of images obtained are then fused following a precise rule before the inverse transformation starts.
The IHS based algorithms have undergone many developments due to the simplicity of their implementation and are often used in the study of QuickBird and IKONOS images. Thus, Tu and et al. (2004) have proposed an extension of the IHS basic algorithm to promote vegetation zones that appear very dark in the green band (G Band) in IKONOS images because of their low power of reflectance. In their algorithm, a fourth component, the near infrared (NIR), is associated with the traditional components of the RGB colorimetric space for the transformation of the MS image in the colorimetric space IHS. The adjustment of the influence of different spectral bands is done thanks to coefficients of spectral adjustment. These are either fixed (Tu et al., 2004), either determined according to the value of parameters such as the Normalized Vegetation Index (NVI) and its variations (Tu et al. 2009;Choi, 2006;El-Mezouar, 2012). The fusion of the images processed can be performed according to the simple rules or much more complex as the use of a genetic algorithm in the Generalized IHS-genetic algorithm (GIHS-GA) (Lari and Yazdi, 2016). These more complex operations permit to optimize the fusion of images.
Hybrid methods have also been proposed in the literature, including those that combine IHS base methods and multire solution analysis. Thus, the Additive Wavelet Proportional Luminance (AWLP) algorithm proposed by Otazu et al. (2005), transforms the MS image in the colorimetric space IHS. The details of the PAN image are then injected in each spectral band in proportionally to a weight which depends on the luminance of the MS original band. In the method described by Lari and Yazdi (2016), the Generalized Laplacian Pyramid with Spectral Distortion Minimization (GLP-SDM), the multiresolution analysis uses a Laplacian pyramid and the injection is carried out according to a weight which depends on the MS image band and the high frequency component of the PAN image (Lari and Yazdi, 2016). The main limitation of the pansharpening methods and in particular those IHS based is that they realize the fusion by taking into account only a part of the information sources. In the classical IHS methods (IHS, FIHS, GIHS), the treatments are only based on an analysis at the level of the pixels taken individually. This analysis is focused upon the intensities of the pixels for the PAN image and the MS image. The different extensions of these methods only take into account the spatial information that only the PAN component contains. This information is usually extracted by a high frequency filtering. Unfortunately, in the images with very high spatial resolution, the increase of the spatial resolution is reflected by a decrease in the spectral resolution. In addition, this type of images presents a large intraclass spectral variability. The consequence of this is that there is a high spectral distortion of the image fusion doubled with spatial distortion, which is characterized by the presence of blurring effects and artefacts. To overcome these limits, some authors (Palsson et al., 2014;He et al., 2014) have recently used the parsimonious analysis of signals which realizes the full extraction of the information contained in an image. The advantage in using the parsimonious analysis to fuse images is that a decimation of the fused image permits to completely rebuild the multispectral image while the panchromatic image can be viewed as a linear combination of the components of the image fusion (Palsson et al., 2014). It leads to a linear system under-determined (the matrix representing the linear system is a rectangular matrix with a more important number of columns than rows) which is resolved by the total variational methods (Rudin et al., 1992). The results are highly superior to the conventional algorithms of image fusions. However, this new approach is limited the one hand by the complexity of implementation of algorithms used (Basis Pursuit, Orthogonal Matching Pursuit, Chambolle-Pock) and on the other hand limited by the need to have a dictionary of functions complete enough. Moreover, the choice of the dictionary is an important factor in getting good performance during the merger of the images. The morphological components analysis (Starck et al., 2005) is a method of parsimonious representation that leans on the morphological diversity of an image and resolves the problems of complexity in algorithm implementations and that of the choice of the dictionary (Xu et al., 2016).
In this study, we propose a new method which belongs to the family of IHS based algorithms and which is based on the morphological components analysis. In the new approach, the MS image is transformed in the IHS colorimetric space before the PAN and MS images are be divided into smooth components and texture parts using the Morphological Components Analysis (MCA). The interest of this decomposition is double. On the one hand it provides descriptors both robust to noise and to the transformations that the images may undergo and on the other hand, it permits to take into account simultaneously the spectral and spatial information contained in each of the images sources. The images of the smooth component are fused with the help of the Laplacian pyramid and bivariate Hahn polynomials are applied on the texture component in order to provide an efficient description of the texture. The image fusion result is obtained by linear combination of the images of the smooth components and texture components from the different treatments.
The rest of the paper is as follows: In section 2 we describe the analysis in morphological components of our approach. Section 3 presents the fusion method proposed. Section 4 is devoted to the phase of experimentation and the presentation of the results obtained. Finally, we terminate this paper by a conclusion.

Decomposition of Images by Morphological Components Analysis (MCA)
An efficient analysis of very high spatial resolution remote sensing images must take into account both the spectral content and the spatial content of the images. Like the conventional approaches commonly encountered in the literature28 and especially in the analysis of the images with very high spatial resolution, we propose to decompose the Panchromatic image (PAN) and Multispectral image (MS) into Smooth component (STRUCT) and Textural component (TEXT). The Morphological Components Analysis (MCA), proposed by Starck and Donoho (2005), offers a parsimonious representation based on the use of dictionaries. It allows to realize the decomposition of PAN and MS images in their components SRUCT and TEXT. According to the formalism of MCA, any image X possibly affected by noise, can be written as a linear combination of its K different morphological components different morphological components {x i | i ∈1,2,…, K} : 1,2,..., Each component x i is parsimoniously represented in a dictionary of Φ i basis of wave atoms. x i = Φ i αi, i ∈ 1,2,…,K where αi is a parsimonious vector that corresponds to the projection of the vector x i in the orthogonal basis Φ i . The choice of the dictionary functions is important for, it affects the reparability level between the different morphological components. In effect, the component x i is sparse in one dictionary Φ i but not in another dictionary Φ j (j ≠i).
The research of the components corresponds to the resolution of the following optimization problem: The search of functions x i = Φ i α i i ∈ 1,2,…,K such as: λ is the constant of penalty usually determined during the algorithm by different thresholding techniques. The most used method is that of the mean of max (MoM method) proposed by Candes et al. (2006) In practice, two types of dictionaries are chosen. The first one is the Undecimated Wavelet Transform (UDWT) for the construction of the STRUCT image (smooth part) because it permits to compensate the lack of invariance to the lost translation observed in the case of the classical wavelet transform. Moreover it is parsimonious for image with isotropic characteristics and non-parsimonious for images with anisotropic characteristics. The second one is the Discrete Cosine Transform (DCT) for the construction of the Component TEXT because it permits to focus the relevant information on a small number of coefficients. The following Fig. 1 shows the STRUCT and TEXT components of a multispectral CEBRS image with 19.5 m spatial resolution and captured in B3 band (0.63 to 0.69 µm), B2 band (0.52 to 0.59 µm) and B1 band (0.45 to 0.52 µm) .

Related Work
The aim of pan sharpening is to combine Pan image and MS image, to produce MS imagery with a higher spatial resolution by using suitable algorithm. The classical pan sharpening algorithms are pixel level fusion techniques. In the case of very high resolution image like QuickBird or IKONOS image, their efficient analysis needs to pay attention simultaneously to their spectral characteristics and their textural characteristics (Chen et al., 2012). This assumption is the main idea that guides the novel approach.

Fusion Rule for Smooth Components
According the general principle of the image fusion techniques, one tries to inject into the MS image the information coming from the high spatial resolution PAN image, without degrading the initial information contained in the MS. The approach used to fuse the smooth components of the two types of images is adapted to our problem by the algorithm proposed in 2012 by Saleem et al. (2012) This algorithm constructs a fused image while improving the contrast. This permits to better distinguish the objects present in the image at the initial level. In this algorithm the fused image is carried out using a multiresolution analysis based on the Laplacian and Gaussian pyramids. The weight matrices relating to each level of resolution are estimated by Gaussian pyramids decomposition of an original matrix. They are defined on the basis of the calculation of metrics upon the contrast and the brightness of the original image.
In our method, we propose to determine the own metrics at each resolution. This leads to some multi resolution metrics that take into account the specific characteristics of the image obtained at a given level of decomposition in the Laplacian pyramids. For the implementation of our algorithm, we determine the matrices of contrast and luminance relating to a given level of decomposition. To build contrast matrix, we use the gradient magnitude. The magnitude of the gradient of an image is strongly correlated to its contrast. Thus, the coefficients of our contrast matrix C are the magnitude of gradient C x,y determined in each pixel (x, y) of a given image I: The improvement of the contrast is reflected by a deterioration of the brightness. To reduce this effect, we determine the coefficients L x,y of the luminance matrix L. These coefficients will be used for weighting the degradation of the brightness of the image I. The coefficients L x,y of the luminance matrix will be determined using a Gaussian kernel. They are estimated as follow: With m I the average and σ I , the variance of the image I. Subsequently, each pixel (x, y) of an image at the resolution d, is affected by the weight , d x y P defined by: This weight reflects the level of contribution of each pixel of images sources in the fused image. And the parameters α and β are true constants. The weight used afterwards will be the normalized one, such as the sum of the components

Fusion Rule for Texture Components
The fusion of the texture components from PAN and MS images is based on a calculation of the texture descriptors in each image. These descriptors are determined through discrete orthogonal moments that are obtained from bivariate Hahn polynomials. These polynomials recently introduced by Wu and Yan (2014) characterize the local variations of an image and the moments deriving from it, are robust to noises. In addition, these moments present a numerical stability, what permits to construct algorithms numerically stable (Wu and Yan, 2014 The bivariate Hahn polynomials are thus written as a tensor product of two Hahn polynomial with one variable x ∈[0; N-1] of order N. The unidimensional Hahn polynomial verify the following recurrence: ( 2)ˆ( , , ), 2,3,..., 1 ( ) x n n N n n n N n n B n n n N n η γ η γ η γ η γ η γ η γ η γ η γ η − + + + + + = + − + + + + + − + + + + + + + + = + + + + + + + − The term (α) x represents the symbol of Pochammer. Constants η and γ influence the form and the symmetry of the Hahn polynomial. The details of the calculation of these polynomials are available in Wu and Yan (2014).
The 25 texture's descriptors built in each pixel represent the correlation between the frequency variations of the region of the considered image and the Hahn basis polynomials. They permit to determine the signature vector of the texture at this pixel. This vector has s = 2n-2 components, with n, as the size of the neighbourhood. This vector is based on the amplitudes of the moments MH pq (.,.) of order s = p+q, that contain the whole information texture. The properties of discrimination of this vector in the classification of textures have been studied by Marcos and Cristobal (2013) for the discrete Tchebychev moments, which are a particular case of the Hahn moments. For n = 5, the signature vector corresponding to each pixel (x, y) has 8 components and it is given by: In a second time, we apply the fusion rule according to the value of the similarity matrix S x,y .
If S x,y ≥ θ then the pixel value (x, y) of the image fusion is given by: The weight being respectively determined by: where, ||.|| 1 represents the standard L 1 norm. If S xy < θ then the value of the pixel (x, y) of the image fusion is obtained by: The threshold θ is estimated by the calculation of the average of the similarities on the vicinity of the pixel (x, y).

Flowchart of Our Method
To take care all information contained in Pan image with high spatial resolution and MS image with high spectral resolution, let us consider the MS image noted I MS represented in true colors (R, G and B) and the high spatial resolution PAN image noted I PAN . Our algorithm evolves according to the following steps: Step 1: Mapping of the MS image with the PAN image by resizing the MS image to the size of the PAN image using the Iterative Algorithm Curvature Based Interpolation (ICBI) proposed by Giachetti and Asuni (2008). The choice of this algorithm is explained by the fact that it eliminates the blur effect introduced by conventional methods while preserving the spectral information (El-Mezouar, 2012 Step 2: We use morphological components analysis to extract smooth part and textural parts of each image. We get on the one hand, the images Step 3: We fuse each morphological component to obtain the grey level fused image. We obtain the smooth and texture fused images, respectively referred to as Step 4: We estimate the difference image between grey level fused image and the intensity component (I) obtained in Step 1 Step 5: The color fused image is reconstruct by linear combination of MS image and difference image the colored image fusion is creating by the help of the following relations which gives the color image fusion matrix. Some coefficients will be injected in accordance with the approach proposed by Tu et al. (2009) and El-Mezouar (2012) in order to adjust our algorithm to the spectral response of the sensors like Quick Bird and IKONOS. So we have: The experiments carried out by El-Mezouar (2012) show that the values of the coefficients ω R = 0.6, ω V = 0.82, ω B = 0.37 give the best results.
The Fig. 2 in the previous page shows the flowchart of the novel method.

Experimentations and Results
The qualitative evaluation of the results of our experiments was conducted by calculating different parameters. These are: (1) The Q4 parameter (Universal Index Quality Image) which combines in a single parameter the correlation between the fused image and the MS image, the spectral distortion between the two images and the contrast distortion, (2) the QNR which assesses the quality of an algorithm of fusion without reference image, (3) the spatial distortion D s , (4) the spectral distortion Dλ, (5) the Bias, (6) the relative global-dimensional error (ERGAS), (7) the spectral Correlation Coefficient (CC), (8) its spatial version (sCC) and (9) the absolute mean Squared Error (RASE). The exact calculation of these different parameters is available in (Lari and Yazdi, 2016). Some of them such as D s and sCC measures the ability of the algorithm to preserve the spatial information contained in the panchromatic image while CC and D λ allow to evaluate the quality of the conservation of the Spectral information contained in the multispectral image. Bias measures the difference between the mean of original image and fused image whereas RASE and ERGAS are global measure index. RASE expressed in percentage characterizes the average performance of the method of image fusion in the spectral bands considered. ERGAS is a global quality index sensitive to mean shifting and dynamic range change. Q4 index measures the local correlation, luminance and contrast between two images whereas QNR index expresses the combination of spectral and spatial distortion.

Fig. 2. Flowchart of the novel method
The smaller the values of the bias, ERGAS and RASE, are, the more effective the algorithm is. In the other direction, the larger the QNR and the Q4 values are, the more the algorithm gives good results. As for the coefficients of correlations, some values close to 1 indicate that the spectral and space information are really quite present in the merged image.
In practice, the parameters η and γ used in the fusion Rule of the texture components, are chosen equal (Wu and Yan, 2014). For our applications, we choose η = γ = 25. The parameters α and β were made respectively equal to 1.5 and 2.5 in our application for the fusion Rule of the smooth components. Our first experiments have concerned the fusion of: • A panchromatic image acquired by the IKONOS satellite with 1 m spatial resolution and 1024×1024 pixels (Fig. 3a) and a multispectral image of 256×256 pixels size provided by the satellite LANDSAT ETM+ (Fig. 3.b) • Two images acquired by the Landsat satellite ETM+ ( Fig. 3c and 3d). The PAN image size is 512×512 pixels and the MS image is of 256×256 pixels size and is represented in false colors In all our experiments, the MS images have been resized to the size of the PAN image using the Iterative Curvature Based Interpolation (ICBI) algorithm proposed by Giachetti and Asuni (2008).
The proposed method has been compared to 3 other methods of images fusion classified among the most efficient (Alparone et al., 2008) in HIS based family algorithms. These are the GIHS algorithms-GA, AWLP and GLP-SDM which have been previously described. Figure 4 shows the results of the image fusion tests and Table 1 contains the values obtained through the assessment of the Q4, D S , D λ and QNR parameters. We can notice that the images obtained with the AWLP show an excellent distribution of colors ( Fig. 4a and 4e) even if they are slightly blurred. This reflects the ability of this algorithm to retain the spectral information. This image is however somewhat slightly blurred. Fig. 4b and 4f indicate that the GLP-SDM really preserves the spectral characteristics of the MS image with a mean level of spatial details conservation. As for the results of the image fusion by GIHS-GA ( Fig. 4c and  4g), the details of spatial and spectral features are quite well drawn but we notice a degradation of colors distribution, which reflects a strong spectral distortion as Lari has noticed (Lari and Yazdi, 2016). The result obtained by our method (Fig. 4d and 4h) shows the presence of a slight effect of blurriness. In contrast, the spectral and spatial details are relatively well preserved.
In addition, the dynamics of the image and therefore the contrast are the best. The details and the forms of our images are much better designed.
The qualitative analysis as shown in Table 1 confirms the visual impression. The results obtained by our algorithm are better as regards to the spectral distortion D λ, D s , QNR and Q4. This indicates that the spectral information contained in the multispectral image is well preserved in the image fusion of our method. The values taken by the Q4 parameter in our case are substantially equal to those of the AWLP method even though this last one presents the best results. As for the spatial distortion, our method is effective on the first type of images and less good on the Land sat images where it nevertheless supplants the GIHS-GA which is the one of the reference algorithms on the IHS family methods.
The effectiveness of our algorithm was also assessed on the two other types of images. The first image shows a rural region acquired by the Quick Bird satellite. This is an image scene on the Kokilai Lagoon, a Marine Protected Area in Sri Lanka, taken by the Quick Bird satellite sensor on April 2005(El-Mezouar, 2012. The image size is approximately 2600×3200 pixels. These images is much more complex with different classes of soil occupation. Fig. 5a and 5b show the panchromatic image and multispectral images. The other images ( Fig. 6a and 6b) show an urban area provided by the IKONOS satellite. This image shows a mixed area combined vegetation areas and urban area. The image size is 256×256 pixel. Our algorithm can be classified among the IHS based methods family and has been here compared to 4 conventional algorithms of this family: The Classic IHS (SA-IHS) and the algorithms proposed by TU (2009;Choi, 2006;El-Mezouar, 2012).
The different results of images fusion are presented in Fig. 7a and 7e for the images provided by the Quick Bird satellite and Fig. 8a and 8e for those from the IKONOS satellite.
The visual analysis of Fig. 7 shows us that the image fused by SA-IHS (Fig. 7c) presents a high spectral distortion compared to the other methods. Our method (Fig. 7g) really preserves the spectral information while improving the contrast of the original multispectral image. The results of Fig. 8 confirm the previous results. The image obtained by our method really retains the spectral information and the contrast is strengthened in it.
The dynamics of the image is more highlighted in our approach compared to the other algorithms ( Fig. 8a and  8e). The algorithm of Choi (Fig. 8c) provides an image fusion with a green excessive hue.   Table 2. They confirm the results of the visual analysis previously observed. The analysis of Table  2 shows that our method is effective at the level of all the quality parameters with the exception of the correlation coefficient for which Choi and El-Mezouar methods give best results, mainly in the Red spectral band.
In table 3, our method remains effective at the level of the quality parameters used in the different spectral bands with the exception of the spectral correlation coefficient. The algorithm SA-IHS gives best results in the red spectral band for the correlation coefficient and surpasses a little bit our algorithm. The algorithm of El Mezouar obtains the best Q4. Moreover, the first results in Bias, Variance and Standard Deviation show that our algorithm does not favour any spectral band in the multi spectral image (IKONOS MS image) in contrast to the algorithms such as those of Choi, TU and El-mezour algorithms which generally favour the Green spectral band by degrading the spectral information contained in the Red and Blue bands.

Conclusion
We have proposed in this study a new method of image fusion for panchromatic and multispectral images. The diversity of the sensors makes it necessary to the development of such methods in order to be able to have a global vision and a synthetic view of the information contained in various data sources. The proposed approach which is a hybrid approach IHS-multiresolution, provides encouraging results, when we see the first experiments. After transforming transformation the original multi spectral image in IHS color space, this novel approach uses morphological component analysis to separate PAN and MS images in smooth part and texture part and fused each morphological parts of them to obtain the fused image. To evaluate the accuracy of our method, we compare performance of the novel approach to the most useful HIS based family algorithms like AWLP, GIHS-GA, GLP-SDM, Tu and Choi algorithms SAIHS and EL-Mezouar algorithms. We use qualitative assessments that allow to predict quality of the algorithm to preserve spatial information containing in panchromatic image and spectral information containing in multi spectral image. For all experiments we made in different type of remote sensing images like QuickBird PAN and MS Image, IKONOS and LANDSAT images. The proposed method gives the best fused images in terms of CC, Q, bias, RASE and ERGAS. The results obtained show a better conservation of the spectral information characterized by a high spectral Correlation Coefficient (CC) which translates a loss of 15% for spectral information compare to those of the conventional methods which revolve around 25%. As for the degradation of spatial information (sCC) it is of the order of 17% in contrast to conventional algorithms which oscillate around 21%. The contrast of the fused image is improved and the MS image is virtually retained. This fact is shown by the smallest value of Q4 (> 70%), QNR (>80%) and ERGAS equals to 3.88 for IKONOS image and 3.65 for QuickBird image compared to 8.70 for IKONOS image and 6.97 for QuickBird for conventional IHS algorithms in all of our experiments.