Investigation on image fusion of remotely sensed images with substantially different spectral properties a

Vladimir Buntilov*a, Timo Bretschneiderb Sch. of Comp. Eng, N4-02a-32, Nanyang Avenue, Nanyang Tech. Univ., Singapore 639798 b EADS Innovation Works Singapore, 41 Science Park Road #01-30, Singapore 117610 ABSTRACT

This paper presents the results of the investigation on image fusion applied to data with significantly different spectral properties. The experiments with simulated and real data show that conventional fusion methods heavily distort geometrical shapes of features in case of fusion of data with different spectral characteristics. As a step forward to solve the problem, a novel fusion approach is proposed. The method uses the relationship between the Wavelet Transform Modulus Maxima (WTMM) corresponding to the same feature in the fusing images. The experiments showed advantages of the proposed method in areas of substantially dissimilar spectral characteristics of merging data. The found limitations of the method are discussed as well.

Keywords: image fusion, pan-sharpening, quality evaluation, wavelet transform modulus maxima, WTMM.

1. INTRODUCTION Due to physical and operational constraints of remote sensing instruments the acquired data might either exhibit high spatial characteristics or maintain a great fidelity of spectral properties of the observed objects. Image fusion techniques aim at providing products, which possess high spatial and spectral accuracy simultaneously. Usually, the highly resolved spatial information from the panchromatic band is merged with the spatially lower resolved multispectral data. The main requirement on such pan-sharpening algorithms is to preserve the original spectral properties of observed objects in the fused product [1], i.e. sharpening should not be performed at the expense of accuracy of radiometric information. Although the necessity for spectral consistency of the fused product has long been recognised by the scientific community [1-7], there are few methods, which are explicitly designed to meet this requirement [8, 9]. If the fusing data have similar spectral characteristics, then simple alternation of the input images, such as histogram matching, help to avoid the problem of colour distortions in the fused product. However, simple modifications are insufficient if a certain sensed object results in significantly different intensity responses in the panchromatic and multispectral channels. It is a well-known fact that in this case the fused output might contain colour distortions [10]. Merging of such type of data with low spectral similarity is still a relatively unexplored issue in image fusion. In this paper the investigation on image fusion applied to data with significantly different spectral properties is presented. The necessity of placing high emphasis on the underscored problem is demonstrated with the use of simulated and real data. The first set of experiments employed artificial one-dimensional profiles with different spatial resolutions and spectral responses of sensors. Manually created data helped to test the fusion methods on simple geometrical objects, which behaviour after the resolution enhancement is well-known or obvious. Moreover, the use of artificial profiles allowed to have the “ideal” profile for the sharpened multispectral data, which greatly simplifies the quality evaluation. Real-world scenarios were tested using SPOT-4 imagery. The data were merged by common fusion methods and several well-known quality measures were calculated to evaluate the degree of introduced distortions. The experiments indicated that conventional fusion methods heavily distort geometrical shapes of features in case of fusion of data with different spectral characteristics. As a step forward to solve the problem, an approach called Content Separation Fusion (CSF) is proposed. The method relies on the separation of spectral and spatial information in the detail coefficients of the wavelet transforms of the ggggggggggggggggggggfffffffffff *[email protected]; phone +66 875 066837 ; fax +66 2 2432437

original images with the help of the Wavelet Transform Modulus Maxima (WTMM) technique. The relationship between the WTMM corresponding to the same feature in the corresponding input images is used to construct the wavelet coefficients of the fused product. Previously the CSF approach was proposed for fusion of general-type optical remotely sensed data [11]. This paper is focused on the application of CSF for fusion of data with substantially different spectral characteristics. The proposed approach was compared with common merging techniques for both artificial and real data. The experiments showed that CSF is able to sharpen the data with substantially dissimilar spectral characteristics. However, the tests on real imagery reveal certain limitations of the proposed technique, which have to be eliminated in order to provide high-quality imagery. The remaining part of this paper is organised as follows: Sec. 2 presents the experiments on spectrally dissimilar data fusion with the simulated profiles of different types of edges. Real-world scenario analogues of fusion are discussed in Sec. 3. The proposed fusion method is introduced in Sec. 4. Finally, Sec. 5 concludes the paper.

2. EXPERIMENTS ON SIMULATED EDGE PROFILES The first set of experiments was conducted using artificial one-dimensional profiles of different geometrical features, namely, step- and roof-edges. All procedures in the experiments were designed in the way to simulate the real-world processes which take place in remote sensing of objects by panchromatic and multispectral scanners. The purpose of using artificially created data was two-fold. Firstly, manually created data allow the availability of the “ideal” signal, i.e. the signal, which the multispectral scanner would produce if it had the sharpness as high as the panchromatic scanner. This “ideal” signal helps to analyse the problems of the fused profiles and significantly simplifies the calculation of quality indices. The second reason for employing synthesised profiles is the possibility to construct geometrical features with simple pre-defined types, such as “pure” step- or roof-edges. The behaviour of such objects under varying parameters (such as spatial resolution, spectral properties etc.) is known and obvious. This helps to focus on the analysis of the introduced distortions and assessment of the actual fusion method. 2.1 Profiles creation Let pan, xs and xs_ideal denote the one-dimensional profiles, which describe the signals obtained by the imaginary high spatial resolution (panchromatic) scanner, lower resolution (multispectral) scanner and the “ideal” high-resolution multispectral scanner, respectively. In the current experiments, each of these profiles has the length of 1000 elements. The emulation of the sensing process of the scanners with different spatial resolutions was performed in the following way. Firstly, the profiles, which represent the real scene “analogous” data before the actual sensing, were constructed. Since the data should be in discrete format for software processing, such profiles are emulated by a deliberately excessive sampling rate: e.g. for each element in pan there are 100 elements in the corresponding “real scene” profile. Thus, the “real scene” profiles, named as s_pan and s_xs, have the length of n=100*1000=106 elements each. A chosen geometrical feature is emulated in each profile at the same position by assigning certain values to the corresponding elements. For example, to emulate a step-edge feature with significantly different spectral responses in panchromatic and multispectral channels, the following Matlab-pseudo code was used: s_pan(1:n/2)=140; s_pan(n/2+1:end)=250; s_xs(1:n/2)=120; s_xs(n/2+1:end)=185;

(1)

The presence of the step-edge produces an increase of 110 values (from 140 to 250) in the emulated panchromatic profile, while for the multispectral profile the increase constitutes 65 values (from 120 to 185). The second step was to emulate the sensing process of the scanners with different spatial resolutions. For this, the profiles s_pan and s_xs were divided into non-overlapped segments with a certain length and the mean of the values inside the segments was calculated, which approximates the sensing and sampling processes of scanners. The actual length of the segments depends on the emulated spatial resolution of the scanners. In the current experiments, the resolution of the imaginary panchromatic instrument was taken as four times higher than the spatial resolution of the multispectral scanner. Thus, to generate the pan profile, the signal s_pan was processed with the segment's length set to 100 elements. A similar procedure was applied to s_xs with the segments length equal to 4*100=400 elements to generate xs. The ideal multispectral profile xs_ideal, which spatial resolution is as high as of pan, was generated by

averaging s_xs within the segments of 100 elements. Since the length of the profiles has to be identical for the fusion, the signal xs was resampled to the size of pan using linear interpolation. The following types of geometrical features were analysed in the experiments. 1.

A step-edge (step_pan_strong), which produces stronger changes in values in pan than in xs due to different spectral responses.

2.

A step-edge (step_xs_strong), which produces stronger changes in values in xs than in pan due to different spectral responses.

3.

A step-edge (step_negative), which produces opposite changes in values in xs compared to pan due to significantly different spectral responses. Such cases, although they might be infrequently in real-world imagery, are still realistic scenarios.

4.

A roof-edge (roof_positive), which produces stronger changes in values in pan than in xs at each of its slopes due to different spectral responses.

5.

A roof-edge (roof_negative), which produces opposite changes in values in xs compared to pan at each of its slopes due to significantly different spectral responses.

6.

A roof-edge (roof_alternate), which produces stronger changes in values in pan than in xs at one of its slopes and opposite changes in values in xs compared to pan at another slope.

The actual values in the first set of profiles (step_pan_strong) are given by (1). The creation of the next two profiles (step_xs_strong and step_negative) is similar and obvious. The profiles of roof-edges were created in a similar way, since a roof-edge can be considered as a combination of two step-edges. The width of the roof-edges was chosen to ensure the existence of an interference effect from the slopes of a roof-edge in a low resolved xs and, as a result, in the fused profiles. Particularly, the width of the step-edges was set to 5 elements in pan and xs_ideal, which corresponds to 500 samples in the s_pan and s_xs profiles. 2.2 Fusion of data The fusion of xs and pan was performed for each type of profile described in Sec. 2.1. In the experiments with the artificial data the emphasis was put on wavelet-based fusion, while other methods were beyond the scope of the current investigation. The purpose of these experiments was to detect and analyse typical problems in fusing data with significantly different spectral properties. Therefore, the simplest well-known wavelet-based merging algorithm was chosen. This ensures that the analysis is focused on typical distortions, rather than on the ability of different fusion rules to improve the results. The profiles xs and pan were firstly decomposed by the wavelet transform up to the third level of decomposition. The approximation of xs provided the approximation of the fused signal fused_WT, the three levels of detail coefficients for fused_WT were taken from the decomposition of pan. Afterwards, the inverse wavelet transform was performed, yielding the actual fused_WT profile. The cubic spline wavelet was used in the experiments with the artificial data, while the details about the decomposition and reconstruction filters can be found in [12]. 2.3 Analysing the results The fusion results are represented in Fig. 1. For the sake of simplicity only the results for step_xs_strong, roof_negative and roof_alternate are shown, since it was found that they sufficiently represent the distortions detected in the conducted experiments. The quantitative analysis of the distortions was performed for each geometrical feature by comparison of fused_WT and xs_ideal profiles with the help of several quality indices, i.e:


Root mean squared error (rmse)




Root mean squared error normalised to the mean value of the ideal signal (rmse_mean)




Correlation(corr)




Quality index by Bovik (q_Bovik) [13]

0

190

-10

d_pan

-20

d_xs

-30

d_xs_ideal

170

pan

160

xs xs_ideal

150

fused_WT

140

Magnitude

Magnitude

180

-40 -50 -60

fused_CSF -70

130

-80

120

-90

110

-100 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

Element number

Element number

(a)

(b) 150

260

pan 240

100

d_xs

fused_WT

200

fused_CSF 180 160

Magnitude

xs_ideal

220

Magnitude

d_pan

xs

d_xs_ideal

50

0

-50 140

-100

120 100

-150 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

Element number

Element number

(c)

(d) 100

260

d_pan

pan 240

xs

d_xs_ideal

fused_WT

200

fused_CSF 180 160 140

Magnitude

Magnitude

d_xs

50

xs_ideal

220

0

-50

-100

120

-150

100 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

Element number

Element number

(e)

(f)

Fig. 1. Fusion results for various edge profiles. (a) step-edge step_xs_strong, (b) the wavelet transforms corresponding to the profiles in (a), (c) roof-edge roof_negative, (d) the wavelet transforms corresponding to the profiles in (c), (e) roof-edge roof_alternate, (f) the wavelet transforms corresponding to the profiles in (e).

All quality indices were calculated for the extracted regions of profiles, in which the corresponding feature shows a strong presence. This ensures that the measures produce objective results and reduces their dependency on the length of the profiles. For example, taking into account the numbering as shown in x-axis in Fig. 1, the quantitative analysis was performed for the elements 10 to 35 for the step-edges and 5 to 40 for the roof-edges. Table I presents the quality indices calculated for all geometrical features employed in the experiments.

Fusion results for the “step_xs_strong” case are shown in Fig. 1a, while the corresponding wavelet transform details for the second level of decomposition are depicted in Fig. 1b. As can be seen from Fig. 1b, the wavelet details of pan (d_pan) are much smaller in values than the ideal details d_xs_ideal (compare black and green profiles). The reason for this is that pan has a smaller jump in values for the step-edge than xs as can be seen in Fig. 1a. As a result, fusing d_pan with the approximation signal of xs yields the profile fused_WT, which shape is heavily distorted compared to xs_ideal. The left and right sides of the fused_WT profile extend even further than the original smoothed xs. Thus, in this case the fusion did not enhance the spatial resolution, but rather distorted and smoothed the original feature. The quality indices in Table I confirm the visual analysis, e.g. rmse value between fused_WT and xs_ideal are larger than between the original xs and xs_ideal. Table I. The quantitative results for the fused profiles of edges.

step_pan_strong

step_xs_strong

step_negative

roof_positive

roof_negative

roof_alternate

xs fused_WT fused_CSF xs fused_WT fused_CSF xs fused_WT fused_CSF xs fused_WT fused_CSF xs fused_WT fused_CSF xs fused_WT fused_CSF

rmse 7.13 10.02 0.41 7.13 8.08 0.41 7.13 39.74 0.41 10.75 11.38 2.51 10.77 45.93 2.46 6.96 20.97 6.92

rmse_mean 0.05 0.07 0.00 0.05 0.05 0.00 0.05 0.26 0.00 0.08 0.09 0.02 0.06 0.26 0.01 0.04 0.12 0.04

corr 0.98 0.99 1.00 0.98 0.98 1.00 0.98 0.08 1.00 0.88 0.99 1.00 0.88 -0.83 1.00 0.99 0.90 0.99

q_Bovik 0.97 0.96 1.00 0.97 0.96 1.00 0.97 0.08 1.00 0.83 0.91 0.99 0.83 -0.80 0.99 0.99 0.90 0.99

The next test case, “roof_negative” is shown in Fig. 1c,d. As can be seen from Fig. 1c, pan has an opposite jump in values near each slope of the feature compared to xs. In real-world imagery this may correspond to a road, which surface material has significantly different spectral responses in the panchromatic and multispectral channels, e.g. in the panchromatic image the road is dark compared to the adjacent areas, but in the multispectral image it is brighter than the areas around it (see the example in Fig. 2). The wavelet details d_pan differ very much from the ideal fused details d_xs_ideal as can be seen from Fig. 1d: d_pan values are not only larger in absolute values than d_xs_ideal, but also they have opposite sign to d_xs_ideal within the extent of the feature. Fusing the profiles with the use of the unmodified d_pan heavily distorts the feature, as can be seen from Fig. 1c. At the same time, Table I shows that rms and rms_mean exhibit manifold increase for fused_WT profile compared to the original xs, while corr and q_Bovik changes their sign. Thus, in the case of roof_negative the fusion process results in a severe degradation of the object. Further utilisation of such fused product for some application, e.g. roads detection, might be difficult. The third example “roof_alternate” presented in Fig. 1e shows the case when the left slope of a roof-edge has a positive change in values in both pan and xs, while the right slope has a positive change in xs and a negative one in pan. Similarly to the previous example, the geometry of the feature is completely distorted and all quality indices for fused_WT reflect a significant deterioration compared to the original xs. In general, as can be seen from Table I, for all profiles analysed in the experiment, the quality indices show a significant or huge deterioration of fused_WT profiles compared to the original xs. This proves that “without a proper modification, fusion of data, which exhibit significantly different spectral characteristics, might result not in enhancement but in a severe deterioration of the geometrical quality of objects.

3. EXPERIMENTS ON REAL DATA The purpose of the second set of experiments is to demonstrate that similar to the discussed in the previous Sec. 2 distortions occur for real data. For this, SPOT-4 multispectral and panchromatic images were selected. The scenes were acquired over Jakarta, Indonesia, in June 2004. Firstly, the four multispectral bands (xs1: 0.50–0.59 µm, xs2: 0.61–0.68 µm, xs3: 0.78–0.89 µm, xs4: 1.58–1.75 µm) with a ground-projected pixel size of 20m×20m were resampled using bicubic interpolation to match the size of the more highly-resolved panchromatic (0.61–0.68 µm) image pan with a spatial resolution of 10m×10m. The method of resampling might have an influence on the final results. However, the goal of the current experiment was to evaluate the distortions produced by well-known fusion techniques. Thus the evaluation did not focus on the resampling technique parameter, while one of the most oftenly used extrapolation method for fusion was chosen, i.e. bicubic resampling. Afterwards, the multispectral bands were co-registered with respect to pan with sub-pixel precision. The following fusion methods were employed in the experiments: 1.

Wavelet-based fusion with full details substitution rule (WT) [2]: Each band of the multispectral image was sharpened individually by combining the approximation of the xsi with the details of pan. The same wavelet as described in Sec. 2.2 was used for analysis. Before the decomposition, histogram matching was applied to pan.

2.

Additive wavelet technique using “a-trous” decomposition method (AWT) [4]: The high-frequency wavelet plane, extracted from the histogram-matched pan, was added to each xsi individually.

3.

Gram-Schmidt spectral sharpening (GS), implemented in the ENVI software package [14]: The Gram-Schmidt transformation was applied to the spatially low resolved image set with the consecutive replacing of the first transformation component by the spatially highly resolved panchromatic data.

4.

PCA-based method (PCA): In this method the first principal component pc1 of xs1..4 was replaced by the histogram matched panchromatic band [15].

5.

IHS-based method (IHS) [15]: In this method the imagery was sharpened by replacing its intensity component I by the histogram matched panchromatic band. Since IHS-based methods can only be applied to triplets of data, the four available bands were recombined into a group of three components xs341, which contained the least correlated band xs3.

6.

Combination of PCA- and wavelet-based techniques (PCA-WT) [16]: In this method the first principal component pc1 of the decorrelated multispectral bands was replaced by its sharpened counterpart pc1* To obtain pc1*, pc1 was spatially enhanced by pan using the WT method.

7.

Combination of IHS- and wavelet-based techniques (IHS-WT) [17]: In this method the intensity component was sharpened by pan using the WT method.

The quality evaluation of the fused products was performed for extracted 100×100 pixels sub-images. Using small subimages allows restricting the fusion performance analysis to particular objects with specific characteristics of interest, i.e. each patch contains a certain feature, which exhibits significantly different spectral characteristics in pan and xs. Hereafter, the extracted sub-scenes are denoted by adding an extra subscript in parenthesis to the band’s name, e.g. the first 100×100 pixels sub-image of xs2 is denoted as xs2(1). The fused products for the first patch are presented in Fig. 2. The shown taxi- and runways of Soekarno-Hatta airport have an opposite spectral response in sub-scene xs3(1) of the near infrared band with respect to pan. As can be seen from Fig. 2.c,d,e the transformation-based methods IHS, PCA and GS significantly change the colour of the taxi- and runways. Methods based on the wavelet transform, i.e. WT, IHS-WT and PCA-WT and AWT, create serious spatial distortions near the edges of the ceiled surfaces and fail to sharpen the imagery. It is worth mentioning that the shown example is a real-world analogue of roof_negative case discussed in Sec. 2. Fig. 3 demonstrates the spatial quality of the fused products more vividly, showing extracted profiles corresponding to the weakly correlated xs3. As can be seen from Fig. 3, WT method does not enhance, but distorts the spatial feature. This example is a real-world case of step_negative profile presented in Sec. 2. Another example for fusion of weakly correlated data is displayed in Fig. 4. With the emphasis on the orange parts around the squared object in the centre, it can be observed that the spectral properties of the sharpened bands are greatly distorted by the IHS, PCA, GS and WT methods. IHS-WT and PCA-WT changed the original colours slightly less, but the distortions are still very noticeable and the results are far from satisfactory.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Fig. 2: Examples for fusion of weakly correlated SPOT-4 data. (a) 100×100 pixels sub-scene xs3(1) of the original near infrared band, (b) corresponding sub-scene of the panchromatic band, (c) IHS, (d) PCA, (e) GS, (f) AWT, (g) WT, (h) IHS-WT, (i) PCA-WT, (j) CSF, (k) IHS-CSF, (l) PCA-CSF.

150 p m3 WT CSF

140

Magnitude

130 120 110 100 90 80 70 1

3

5

7

9 11 13 15 17 19 21 23 25 27 29

Pixel number

Fig. 3: Cross-section profiles of step-edges from xs3 SPOT-4 fused product.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Fig. 4: Examples for fusion of weakly-correlated data . (a) 100×100 pixels sub-scene xs341(2), (b) corresponding sub-scene of the panchromatic band, (c) IHS, (d) PCA, (e) GS, (f) AWT, (g) WT, (h) IHS-WT, (i) PCA-WT, (j) CSF, (k) IHSCSF, (l) PCA-CSF.

4. PROPOSED CSF METHOD The experiments in Sec. 2 and Sec. 3 clearly show a need in developing a fusion method, which is able to deal with a problem of merging data with significantly different spectral properties. Conventionally, wavelets are used for low-level content separation: the high-frequency coefficients correspond to the spatial information, while the low-frequency coefficients contain the spectral-dependent content. Therefore, wavelets were chosen as a starting point to overcome the problem of spectrally dissimilar data fusion. However, the carried out in Sec. 2 analysis showed that the traditional wavelet transform is only partially able to split the information into spatial-dependent and spectral-dependent components. Hence, the solution for this problem lies in a further separation of spectral- and spatial-dependent information within the sets of wavelet coefficients. In the proposed novel fusion method, which is named Content Separation Fusion (CSF), “pure” spectral and spatial signatures of visible objects are extracted before the actual fusion. In order to accomplish this task, the wavelet transform modulus maxima technique is employed. Afterwards, the signatures from the panchromatic and multispectral bands are combined and used to reconstruct the actual fused high-frequency wavelet coefficients. The method was firstly introduced by the authors for the merging of remotely sensed data in [11]. Current investigation emphasises the applicability of CSF to one- and two-dimensional data with substantially different spectral characteristics. 4.1 Spatial-dependent information Edges of physical objects on Earth are represented in the sensed data as sharp variations in pixels’ intensities convolved by the imaging system’s inherent point spread function. Hereafter, a one-dimensional spatial feature is defined as an interval of pixels where the data has certain variability. This definition can be straightforwardly extended to the twodimensional case. Thus, wavelets are especially useful for describing spatial features since the detail coefficients tend to have non-zero values where variations occur within the scene. To characterise the spatial-dependent information of features, their position and shape descriptors are established in the following manner: It is natural to define the position of a spatial feature as the location where it locally exhibits the largest degree of variability. In the proposed method, the high-frequency wavelet coefficients are used as an analogue of the gradient. For further details, the work by Mallat et al. provides more theoretical considerations regarding the description of a signal’s variability by its undecimated wavelet transform [18]. In summary, it was shown that the wavelet transform modulus maxima (WM) occur at points of locally extreme variability, provided that the wavelet function is chosen appropriately [18]. Therefore, in this work, the positions of features at a given scale of decomposition are defined as the positions of corresponding wavelet transform modulus maxima. The next step is to represent the shape of a feature from its wavelet transform as the second descriptor of spatialdependent information. Obviously, the actual shape of a feature is encompassed in both the low and high frequency components of the signal. The low-frequencies are automatically taken into account in the approximation part of the wavelet transform and, for the fused image, they are usually taken “as is” from the multispectral band. Thus, hereafter only the high-frequency part of a shape’s description is considered in more detail. Mathematically, an irregularity in a signal can be characterised by its Lipschitz exponent, which reflects the order of the polynomial, approximating the signal within the neighbourhood of the irregularity [18]. The rate of how the WMs, corresponding to the same feature, change relatively from scale to scale can be used to describe the sharpness of a feature [12, 18]. Thus, in the current work this rate is referred to as the high-frequency part of a feature’s shape description. As shown later in Sec. 4.3, the proposed fusion method does not require explicit calculation of the features’ Lipschitz regularity, which significantly simplifies the implementation. In summary, the high-frequency part of spatialdependent information of a sharp variation can be characterised by the locations of the corresponding wavelet transform modulus maxima and their evolution through different scales. 4.2 Spectral-dependent information The digital number (DN) of a pixel is linearly related to the integral over the energy reflected and emitted from the corresponding Earth surface within a given spectral range. By virtually increasing the spatial resolution and keeping the spectral sensitivity untouched, the total energy received by the pixels corresponding to the analysed area remains identical. The total energy sensed by the scanner from the surface of the feature can be represented by the area under the feature’s profile in the image, and is for both the low- as well as highly-resolved case identical. Accordingly, the spectral

information of a feature is sufficiently described by the homogeneous DNs of the pixels before the feature’s spatial occurrence and the difference of DNs after the occurrence. The outlined approach coincides with the wavelet transform splitting the original signal into a low-frequency component, which shows a general trend of DNs at a coarse level, and high-frequency components representing the variability of the data. It has been found that the magnitude of the wavelet transform modulus maxima is directly related to the DNs’ variability within the neighbourhood of a corresponding singularity, provided that the transform is performed using a proper type of wavelet [12]. The width of the wavelet can be altered by changing the decomposition level of the analysis. Consequently, the spectral-dependent content of features can be denoted by their WM values at the analysis level of decomposition. 4.3 Fusion of spatial- and spectral-dependent content In the proposed technique, the individual spatially low-resolved bands are sharpened by the more highly-resolved band. The source with higher confidence in spatial content, i.e. the panchromatic image, is used to extract the positions and shapes of spatial features, while the spectral-dependent information is taken from the multispectral bands. The extracted spectral and spatial information are fused and the result is represented by a set of WM, which are used to construct the actual wavelet details employing an iterative reconstruction algorithm [12]. After the pre-processing of the images, i.e. resampling and alignment, the undecimated wavelet transform is applied using quadratic spline wavelet filters [12]. The deepest level of decomposition, which in this work is called the analysing level, is used to extract the spectral-dependent information of the features. The actual level, i.e. a level at which the WM of xs and pan are compared to determine their spectral relationship, is chosen according to the spatial resolution ratio of the image set. For example, in case of SPOT-4 imagery, the panchromatic image has 10m×10m resolution while the multispectral bands have 20m×20m spatial resolution. The first wavelet decomposition reduces the resolution of pan twice. Hence at decomposition levels deeper than the first, the spatial resolution of pan and xs of SPOT-4 imagery can be considered identical. As a result, the second level of decomposition is used in case of SPOT-4 imagery since deeper levels lose the spatial details. After the wavelet transformation of xs and pan, the wavelet transform modulus maxima of the horizontal (H) and vertical (V) components WM 1{..HN,V }({xs , pan}) are detected. As was shown in Sec. 4.1, the spatial-dependent information of a feature is defined as the positions of WM points corresponding to this feature at different scales and their relative changes from scale to scale, which accounts for the Lipschitz regularity of a variation. Since the spatial information of the fused feature is taken from pan, it is reasonable to transform modulus maxima of pan using a linear model: a

WM 1{..HN,Va −}1 ( f )

s

= {sH,V } ∗ WM 1{..HN,Va −}1 ( pan) ,

(2)

s

where s denotes the detected spatial feature observable in both xs and pan. In this case the spatial information from pan is preserved since this operation does not create new maxima. Moreover, relative changes across scales remain untouched, thus, keeping the feature’s Lipschitz regularity unchanged without explicitly calculating it. The scaling
coefficients are determined for each spatial feature individually at the analysing level Na. The coefficients {sH,V } are {H ,V } {H ,V } chosen in such a way that WM N a ( f ) = WM N a (xs ) , i.e. the spectral-dependent properties of the fused feature are s s the same as of the original one at the analysing level:  H,V (3) = WM N{H ,V } (xs ) WM N{H ,V } ( pan ) . s a

a

S

A few words should be said regarding the construction of WMs if the spectral coefficient for the corresponding feature is absent. This might happen if the algorithm fails to establish the correspondence between the extrema of xs and pan for a given feature at Na. In this case Equations (2) and (3) are not applicable and there is a decision to be taken how to construct the wavelet extrema of the fused image for this feature. Two most obvious solutions are either to take the maxima from WM1.. N a −1 ( pan ) or from WM1.. N a −1 (xs ) without any changes. The first case corresponds to the reconstruction of the feature as it appears in pan, the second one corresponds to the reconstruction of the feature as it appears in xs. However, one of the constraints which are usually imposed on pan-sharpening algorithms is to preserve spectral-dependent characteristics of features, which are contained in the multispectral image. Thus, inserting WM1.. N a −1 ( pan ) may enhance the spatial resolution of the image, but at the same time might degrade the spectral quality of the fused product. Moreover, even solely spatial enhancement is not guaranteed in case the spectral characteristics of a

feature are very different in xs and pan, this issue was discussed in detail with examples in Sec. 2-3.Therefore, it was decided to use the extrema of the multispectral band for the fused image, i.e. the feature in the fused band will appear similar to as it does in the multispectral image. This at least ensures that the spectral properties are kept. There is a drawback of this approach, which is particularly apparent for a human observer: there are some regions in the fused image which look as blurred as in the multispectral image. Such regions bordering with the spatially enhanced areas create an unnatural visual appearance. More details regarding the actual implementation of CSF are given in [11].

5. EVALUATION OF THE PROPOSED CSF METHOD The evaluation of the proposed fusion method was performed for both artificial and real data examples presented in Sec. 2-3. 5.1 CSF performance evaluation for the artificial edge profiles

The experiments with the artificially created edge profiles show a significant superiority of CSF over other tested techniques. In Fig. 1a the profile of step_xs_strong corresponding to CSF fusion is almost indiscernible from the ideal profile. For roof_negative case, the fused by CSF signal is very close to xs_ideal, contrary to the WT fused profile. The shown in Fig. 1e results for roof_alternate represent the case when CSF detects that it is not able to sharpen the data, thus the algorithm reconstructs the original low-resolution data (this behaviour was discussed in Sec. 4.3). The reason for this is that the algorithm detects one extremum in d_xs and fail to find the corresponding extremum in d_pan, since it is located not within the analysed proximity (which is two elements to the left and right). As can be seen, the fused by CSF profile coincides with the original xs, which, in this case provides a better approximation of the ideal signal than the totally distorted profile obtained by WT fusion. The quantitative evaluations presented in Table I show that fused_CSF profiles are closer to xs_ideal than fused_WT for all tested examples. For some cases (e.g. step_negative and roof_negative, the quality indices for fused_CSF show huge improvement compared to fused_WT. For all examples, except roof_alternate, CSF improves the quality of the edges, while for roof_alternate CSF produces the original xs. Contrarily, the analysed WT method degrades the geometry of the objects for all cases. 5.2 CSF performance evaluation for real remotely sensed imagery

Besides the original CSF, the combination methods IHS-CSF and PCA-CSF were tested. These combination methods are similar to IHS-WT and PCA-WT described in Sec. 3, but for IHS-CSF and PCA-CSF the CSF is used instead of WT. The experiments with real imagery reveal both advantages and limitations of the proposed CSF algorithm. As can be seen from Fig. 2j,k,l, the intensity values and geometry of the taxi- and run-ways are preserved and the object became sharper. Among those methods, which more or less preserve spectral characteristics, the visually sharpest results are produced by the methods which incorporate the proposed CSF technique, namely, CSF, IHS-CSF and PCA-CSF. However, the sharpness of the CSF-based products is not as high as of the panchromatic image. From analysing the profiles in Fig. 3 it can be concluded that the CSF-based methods are more suitable than the WTbased counterparts for this case: the latter do not enhance, but distort the spatial feature. The colour examples in Fig. 4j,k,l show similar trends for CSF: the colour preservation of the product is very good, but the spatial resolution is less than of the panchromatic band, although the sharpness increases. It must be noted that the following limitations of the CSF method have to be eliminated before it can be used for realworld applications in the industry. The proposed method has to establish the correspondence between features in pan and xs in order to enhance the spatial information. Thus, if it fails to find this relation, only the spectral and spatial information from xs contribute to the fused product. Therefore, the images sharpened by CSF-based methods might contain problematic areas that totally coincide with the low resolved xs. This unnatural patchy look is unfamiliar to the human observer, especially when a prolonged physical object extends over both problematic and sharpened regions as shown in Fig. 5a,b,c (initial experiments with IKONOS-2 imagery). Another reason for the degradation of CSF-based fusion products is that the degree of accuracy, with which the spatial descriptors of a feature are determined, depends on the presence of other objects in the vicinity. The best accuracy is

achieved in case of an isolated object. However, a strong influence of adjacent features changes the model of WM propagation across scales [18], thus, making the proposed fusion model described by Equations (2) and (3) incapable. Whether a spatial feature can be regarded as isolated or not was discussed in [19]. Therefore, the CSF-based methods show poor performance in areas in which spatial features are located close to each other, e.g. the texture region shown in Fig. 5d,e,f. A few gates of the airport are located close to each other, and their structure is better represented in the more highly resolved panchromatic band of SPOT-4 (Fig. 5e) than in the low resolved multispectral bands (Fig. 5d), which spatial resolutions are not sufficient to capture the fine details. In this case, the CSF technique fails to successfully recover the adjacent objects, which can be seen in Fig. 5f.

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 5: Example of scenes where CSF fails to enhance the imagery. (a) The colour combination of the IKONOS-2 multispectral bands, (b) Panchromatic band corresponding to (a), (c) Upper-right part of the rectangular object in the centre is successfully sharpened by CSF, while the bottom-left part is blurred (taken without changes from multispectral image), (d) False colour combination of SPOT-4 multispectral bands, (e) Panchromatic band corresponding to (d), (f) CSF method distorts the textured object in the centre.

6. CONCLUSION This paper presents the results on the investigation of fusion of data, which exhibit significantly different spectral characteristics. The experiments employed both artificially created edge profiles of various types as well as real remotely sensed imagery of SPOT-4 panchromatic and multispectral instruments. The use of the artificial data allowed to carefully analyse the problems which common fusion methods produce when dealing with signals of substantially different spectral properties. It has been proven that merging the data without a prior proper adjustment will result in significant distortions not only in colours, but also in geometry of the sensed objects. The proposed Content Separation Fusion method (CSF) tries to eliminate the problem of dissimilar data merging. It was assumed that the spatial information in the fused product should be taken from the highly spatially resolved data (e.g. panchromatic band), while, at the same time it should be modified according to the spectral characteristics of the less sharp data (multispectral band). To accomplish this, the Wavelet Transform Modulus Maxima technique (WTMM) was employed to separate spatial and spectral information in the wavelet transform of the signals. The analysis of CSF-based fusion results reveals both its advantages and limitations. The proposed method is able to sharpen the objects with significantly different spectral characteristics in the fusing images. For such type of features CSF often outperforms other tested merging methods. However, the products of CSF are not as sharp as the panchromatic data. In addition, other discussed limitations of the CSF method have to be eliminated before it can be used for real-world applications in the industry.

REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]

[11]

[12] [13] [14] [15]

[16] [17] [18] [19]

L. Wald, T. Ranchin, and M. Mangolini, “Fusion of satellite images of different spatial resolutions: assessing the quality of resulting images,” Photogrammetric Engineering and Remote Sensing 63(6), 691–699 (1997). D. Yocky, “Image merging and data fusion by means of the discrete two-dimensional wavelet transform,” Journal of the Optical Society of America A 12(9), 1834–1841 (1995). A. Garzelli, and F. Nencini, “Interband structure modeling for Pan-sharpening of very high-resolution multispectral images,” Information Fusion 6(3), 213–224 (2005). J. Nunez, X. Otazu, O. Fors et al., “Multiresolution-based image fusion with additive wavelet decomposition,” IEEE Transactions on Geoscience and Remote Sensing 37(3), 1204–1211 (1999). Z. Wang, D. Ziou, C. Armenakis et al., “A comparative analysis of image fusion methods,” IEEE Transactions on Geoscience and Remote Sensing 43(6), 1391–1402 (2005). J. Zhou, D. L. Civco, and J. A. Silander, “A wavelet transform method to merge Landsat TM and SPOT panchromatic data,” International Journal of Remote Sensing 19(4), 743–757 (1998). J. C. Price, “Combining multispectral data of differing spatial resolution,” IEEE transactions on Geoscience and Remote Sensing 37(3), 1199–1203 (1999). H. Aanaes, J. R. Sveinsson, A. A. Nielsen et al., “Model-Based Satellite Image Fusion,” IEEE Transactions on Geoscience and Remote Sensing 46(5), 1336–1346 (2008). J. G. Liu, “Smoothing filter-based intensity modulation: a spectral preserve image fusion technique for improving spatial details,” International Journal of Remote Sensing 21(18), 3461–3472 (2000). F. Laporterie-Dejean, H. de Boissezon, G. Flouzat et al., “Thematic and statistical evaluations of five panchromatic/multispectral fusion methods on simulated PLEIADES-HR images,” Information Fusion 6(3), 193– 212 (2005). V. Buntilov, and T. Bretschneider, “A Content Separation Image Fusion Approach: Towards Conformity between Spectral and Spatial Information,” IEEE Transactions on Geoscience and Remote Sensing 45(10), 3252–3263 (2007). S. Mallat, and S. Zhong, “Characterization of signals from multiscale edges,” IEEE Transactions on Pattern Analysis and Machine Intelligence 14(7), 710-732 (1992). Z. Wang, and A. C. Bovik, “A universal image quality index,” IEEE Signal Processing Letters 9(3), 81–84 (2002). C. A. Laben, and B. V. Brower, [Process for enhancing the spatial resolution of multispectral imagery using pansharpening] Eastman Kodak Company, USA(2000). V. K. Shettigara, “A generalized component substitution technique for spatial enhancement of multispectral images using a higher resolution data set,” Photogrammetric Engineering and Remote Sensing 58(5), 561–567 (1992). J. Li, Y. Zhou, and D. Li, “PCA and wavelet transform for fusing panchromatic and multi-spectral images,” Proceedings of the SPIE 3719, 369–377 (1999). A. S. Kumar, B. Kartikeyan, and K. L. Majumdar, “Band sharpening of IRS-multispectral imagery by cubic spline wavelets,” International Journal of Remote Sensing 21(3), 581–594 (2000). S. Mallat, and W. L. Hwang, “Singularity detection and processing with wavelets,” IEEE Transactions on Information Theory 38(2), 617-643 (1992). L. M. Bruce, “Isolation criteria for the wavelet transform mod-max method,” IEEE Transactions on Circuits and Systems II: Analog and Digital Signal Processing 45(8), 1084–087 (1998).