Visible， infrared， and infrared polarization images individually captured by different sensors present complementary information of the same scene， and they could be combined by the image fusion technology to obtain a new， more accurate， comprehensive， and reliable image description of the scene ^［1］. The fusion methods are varied with image sources， fusion requirements， or purposes ^［2-4］. In general， the fusion methods can be classified into pixel-， feature- and decision-level. Compared with the latter two fusion levels， the first fusion level can maintain the source image data as much as possible， so it plays an important role in most image processing tasks. Major pixel-level image fusion methods can be put into four groups according to their adopted theories ^［5］， namely multi-scale decomposition-based methods， the sparse representation-based methods， methods in other domains， and methods combining different transforms. For the multi-scale decomposition-based methods the decomposition schemes and fusion rules are two aspects that affect fusion quality and efficiency.

For the decomposition schemes， various methods are proposed， like the discrete wavelet transform （DWT） ^［6］， dual-tree complex wavelet transform （DTCWT） ^［7］， stationary wavelet transform （SWT） ^［8］， wavelet packet transform （WPT） ^［9］， non-subsampled contourlet transform （NSCT） ^［10-11］， and non-subsampled shearlet transform （NSST） ^［12］. And many practices proved that NSCT and NSST usually outperform other multi-scale decomposition-based methods in representing 2-D singular signals contained in digital images ^［13］. But the design of multi-directional filter banks for NSCT and NSST is relatively complex and computational time-consuming， which greatly reduces the efficiency of image fusion.

Fusion rules generally include low and high frequency coefficients fusion rules. The AVG-ABS rule is a simple fusion rule which uses the average rule to combine low-frequency coefficients and uses the absolute maximum rule to combine the high-frequency coefficients. The AVG-ABS fusion rule is computed easily and implemented simply， however， it always causes distortions and artifacts ^［14-15］. To overcome these shortcomings and improve the fusion quality， a large number of rules have been proposed ^［15-20］. These rules in Refs.15-20 have achieved satisfactory results， but they have the disadvantage of high computational complexity.

To ensure both the fusion quality and computational efficiency simultaneously， a novel multi-scale decomposition-based fusion method with dual decomposition structures is proposed. Our method is dedicated to improving the image fusion quality and efficiency from the aspect of image decomposition scheme， while for the rule aspect， our method only uses simple AVG-ABS rule. Firstly， inspired by the idea of constructing octaves in SIFT ^［21］ and SURF ^［22］ algorithms， the source images are decomposed into a series of detail and approximation images by multi-scale Gaussian filters to construct the undecimated pyramid structures. The multi-scale Gaussian filters have increasing standard deviation as well as up-scaling size. Secondly， for the approximation images， i.e.， the top layers of the undecimated pyramid structures， multi-scale morphology top- and bottom-hat decompositions ^［23-24］ are used to fully extract bright and dark details of different scales on the background， and then the contrast of the fused layer is improved by the absolute maximum rule. Thirdly， the multi-scale morphological inner- and outer-boundary decompositions are especially constructed based the idea of constructing multi-scale top- and bottom-hat decompositions. For each detail image， these two morphology decompositions are implemented to extract the boundary information. And then the decomposed coefficients are combined by the approach of choosing absolute maximum. At last， the fused image is reconstructed through taking the inverse transforms corresponding to the decompositions mentioned.

The theory and mathematical representation for constructing multiresolution pyramid transform scheme are presented by Ref.25 and extended by Ref.26. A domain of signals is assigned at each level， the analysis operators maps an image to a higher level in the pyramid， while the synthesis operator maps an image to a lower level in the pyramid， i.e. and . The detail signal contains information of x which does not exist in ， where and is a subtraction operator mapping into the set . The decomposition process of an input image f is expressed as Eq.1：

f0→y1,f2→…→y1,y2,…,yj,fj+1→…

where

f0=f∈V0fj+1=ψj↑fj∈Vj+1,yj=fj-ψj↓fj+1∈Yjj≥0

And the reconstruction process through the backward recursion is expressed as Eq.3：

f=f0,fj=ψj↓fj+1+yj,j≥0

Eq.1 and Eq.3 are called the pyramid transform and the inverse pyramid transform respectively.

The scale space of an image can be generated through convolving the image with Gaussian filters， and it has been successfully applied in SIFT ^［21］ to detect key points which are invariant to scales. In Ref.26 the scale space is divided into octaves. For each octave， the initial image is iteratively convolved with Gaussians with increasing standard deviation to generate a set of scale space images （Gaussian images）， and one of the Gaussian images is downsampled to obtain the initial image of the next octave. Then， the Difference of Gaussians （DoG） images are obtained by subtracting adjacent Gaussian images. In SURF ^［22］， in order to omit the down-sampling step， the scale space is obtained by increasing the size of filter.

Inspired by the above algorithms， we have the source image repeatedly convolved with Gaussian filters whose standard deviation and size increase simultaneously to construct undecimated pyramid structure. Then， the DoG images are produced by subtracting adjacent Gaussian images. Accordingly， the transform scheme of such pyramid is given by Eq.4：

f0=f∈V0fj+1=Gj*fj∈Vj+1,yj=fj-fj+1∈Yjj≥0

where

Gj=1/2πσj2∙exp-x2+y2/2σj2

is the Gaussian kernels （filters） with a size of in this paper， and denotes the convolution operation. The parameter is the standard deviation， which is increasing with ， and in this paper . Then the source image f can be decomposed into an approximation image and a set of detail images as shown in scheme 1， and it can also be exactly reconstructed through the following recursion：

f=f0

The four-level decomposition scheme is illustrated in Fig.1.

The multi-scale top-hat transform using structuring elements with up-scaling size can extract the light and dark details at different image scales in image fusion ^［24］. Based the idea of constructing multi-scale top-hat transform， the multi-scale morphological inner-boundary transform is constructed. These two kinds of morphological transforms can be expressed as Eq.4， with the Gaussian kernel being replaced by morphological opening operation and erosion operation ， respectively. For purpose of extracting details of different scales， the scale of structuring element increases with j. The inverse transforms can be expressed as Eq.6.

The multi-scale morphological bottom-hat transform and its inverse are shown as follows

f0=f∈V0fj+1=Closebj(fj)∈Vj+1,yj=fj+1-fj∈Yjj≥0

f=f0,fj=fj+1-yj,j≥0

where the analysis operator is morphological closing operation ， with also increasing with . The morphological outer-boundary transform and its inverse are similar to the bottom-hat transform and its inverse， with being replaced by dilation operation .

The proposed fusion method comprises three processes that are multi-scale decomposition， fusion， and reconstruction.

The K-level decomposition of a given source image by the scheme （4） has the form

f→y1,y2,…,yk,…,yK,f(K+1)

where represents the detail image at level and denotes the approximation image of this multi-scale structure.

is a coarse representation of and usually inherit a few bright and dark details， thus the multi-scale top- and bottom-hat decompositions are used to extract bright objects on a dark background and dark objects on a bright background of different scales， respectively. Henceforth， can be decomposed by schemes mentioned in subsection 1.3 as

fK+1→yt(1),yt(2),⋯,ytl,⋯,yt(M),ft(M+1)yb(1),yb(2),⋯,ybl,⋯,yb(M),fb(M+1)

where and represent the detail images at level l obtained by the top- and bottom-hat decomposition process， respectively. And and denote the approximation images of the multi-scale top- and bottom-hat structure， respectively. Figure 2 is given as an example of three-level top- and bottom-hat decompositions.

The detail image in scheme 9 comprises various details like edges and lines， thus the multi-scale inner- and outer-boundary transforms mentioned in subsection 1.3 are used to extract inner-boundary as well as outer-boundary information of different scales. Hence， can be decomposed as

yk→yi(k,1),yi(k,2),⋯,yik,l,⋯,yi(k,Nk),fi(k,Nk+1)yo(k,1),yo(k,2),⋯,yok,l,⋯,yo(k,Nk),fo(k,Nk+1)

where and represent the detail images at level l of that are obtained by the inner- and outer-boundary decomposition process， respectively. and are the approximation images of at the highest level of the multi-scale inner- and outer-boundary structure， respectively. Figure 3 gives an example of three-level inner- and outer-boundary decompositions.

In this paper， the composite approximation coefficients of the approximation image in the multi-scale top- and bottom-hat structures take the average of the approximation of the sources. For the composite detail coefficients of the detail images， the absolute maximum selection rule is used.

The vector coordinate is used here to denote the location of an image. For instance， represents the detail coefficient for the multi-scale top-hat structure at location within level l of source image A. And the notation will be used to denote an image， e. g.， refers to the detail image.

The arbitrary fused detail coefficient and the fused approximation coefficient of the multi-scale top-hat structure are obtained through

ytln|F=maxytln|A,ytln|BftM+1n|F=αtftM+1n|A+βtftM+1n|B

The weights and take 0.5， which preserves the mean intensity of the two source images. Likewise， and of the multi-scale bottom hat structure are obtained through

ybln|F=max{ybln|A,ybln|B}fb(M+1)n|F=αbfb(M+1)n|A+βbfb(M+1)n|B

with .

The selective rule in Eq.12 means that we choose the brighter ones in the bright details， and the selective rule in Eq.13 means that we choose the darker ones in the dark details. In this way， the bright and dark details of different scales can be fully extracted and hence the contrast at each layer can be improved.

For an arbitrary fused detail coefficient of the multi-scale inner-boundary structures， we only use the absolute maximum selection rule：

yi(k,l)n|F=yi(k,l)n|Aifyi(k,l)n|A>yi(k,l)n|Byi(k,l)n|Botherwise

So is the fused approximation coefficient . In such way， the boundary information such as edges and lines of different scales can be well preserved. Likewise， arbitrary and of the multi-scale outer-boundary structures are also obtained by the absolute maximum selection rule.

According to Eqs.6 and 8， the reconstruction of the approximation image can be obtained through the multi-scale top- and bottom-hat inverse transforms as

fK+1∙|F=γM+1ft(M+1)∙|F+∑l=0Mγlytl∙|F/2+κM+1fb(M+1)∙|F-∑l=0Mκlybl∙|F/2

which means both bright and dark information are of equal importance to the source image. In addition， we attach equal importance to the features of different scale levels， thus the weights in Eq.15 are set to be .

Similarly， inner- and outer-boundary information are considered to be equally important to the source image， and so are the features of different scale levels. Thus， according to Eqs.6 and 8， the reconstruction of an arbitrary detail image through the multi-scale inner- and outer-boundary inverse transforms can be obtained as

yk∙|F=fik,Nk+1∙|F+∑l=1Nkyik,l∙|F/2+fo(k,Nk+1)∙|F-∑l=1Nkyok,l∙|F/2

At last， the fused image can be reconstructed by

f∙|F=f0∙|F=fK+1∙|F+∑k=1Kyk∙|F

In order to validate the performance of the proposed method， experiments are conducted on two categories of source images including ten pairs of infrared-visible images （Fig.4（a）） and eight pairs of infrared intensity-polarization images （Fig.4（b））. The two source images in each pair are pre-registered and the size of each image is set to 256×256 pixels. The experiments in this paper are programmed by Matlab 2016b and run on an Intel（R） Core（TM） i5-6500 CPU @ 3.20GHz Desktop with 16.0 GB RAM.

Various pixel-level multi-scale decomposition-based methods including DWT， DTCWT， SWT， WPT， NSCT， and NSST are compared with the proposed method. All the compared methods adopt the simple AVG-ABS rule. According to Ref.13， most of the methods mentioned above perform well when the decomposition levels for them are set to 3. Thus， for purpose of making reliable and persuasive comparisons， the decomposition levels for the methods mentioned above are all set to 3. And to make each method achieve a good performance the other parameters are also suggested by Ref. 13， some of which are listed in Table 1.

For NSST， the size of the local support of shear filter at each level are selected as 8， 16， 32. As for the proposed method， the parameters and k for the multi-scale Gaussian filtering process in Eq. 5 are selected experimentally. In this experiment， the source images are decomposed by 3-layer multi-scale Gaussian decomposition， and different fused images are obtained by changing the parameters and k. During the fusion process， the AVG-ABS rule is also adopted. When and k are in certain value， every fusion image will be evaluated by seven objective assessment metrics （mentioned in subsection 3.2）. For each metric， its mean value is obtained by averaging the evaluation results of the fusion images. Then， seven mean values are summed to get the sum values of objective metrics. Figure 5 gives three surface plots which show variations of the sum of the seven metrics with and k. As shown in Fig.5， the optimal values of and k for the four kinds of images are obtained. The structuring elements in the multi-scale inner- and outer-boundary decompositions are selected as square， and in the multi-scale top- and bottom-hat decompositions they are chosen to be disk. and k in Eq. 5 and the parameters K， M， N₁， N₂， and N₃ in schemes 9， 10， and 11 are set as shown in Table 2 to make the proposed method achieve a good performance.

Seven representative metrics， i.e.， Q₀^［27］， Q_E^［28］， Q^AB/F［29］， information entropy （IE） ^［30］， mutual information （MI） ^［31］， Tamura contrast （TC） ^［32］， and visual information fidelity （VIF） ^［33］ are employed to evaluate the proposed method comprehensively. The variable in TC is chosen to be 4.

In this section， the subjective assessment of the fusion methods is done by comparing the visual results obtained from the above and proposed methods. One sample pair in each type of source images are selected for visual comparison as shown in Figs.6 and 7.

In Fig.6， both the DWT and WPT methods distort the edges of the roof， which was shown clearly in magnified squares. The DTCWT， SWT， NSCT and NSST methods produce artificial edges in the sky around the roof， while the result obtained by the proposed method is free from such artifacts or brightness distortions. In addition， the walls and the clouds in the sky in Fig.6（i） are brighter those in Fig.6（g） and （h）， which means that the fused image of the proposed method has better contrast.

The edges of the car are distorted heavily in Fig.7 （f）， and slightly distorted in Figs.7（c-e） which is shown more clearly in the corresponding regions in magnified square. And Figs.7（c-h） show some artifacts around the edges of the car. However， in Figs.7（i） there are no distortions or certain artifacts. In addition， the car in magnified square of Fig.7（i） is the darker than those in Figs.7（h） and （i）， which demonstrate that the proposed method has better contrast.

The above experiments confirm that the proposed method performs better in visual effect for the two categories of source images. Although adopting the simple AVG-ABS rule， the proposed method does not generate certain artifacts or distortions and simultaneously preserves the detail information of source images as much as possible.

The objective assessment of the seven multi-scale decomposition-based methods are shown in Tables 3. For the infrared-visible images， the proposed method performs the best on all the seven metrics. For the infrared intensity-polarization images， the proposed method performs the best on the other five metrics except Q₀ and Q_E on which it performs the second best. It can also be obtained from Tables 3 that compared with the seven methods， the proposed method always has the best assessment on metrics Q^AB/F， IE， MI， TC， and VIF. It means that the proposed method can transfer the original information of source image including the edges and brightness details to the fused image sufficiently， and improve the contrast of the fused image.

To verify the efficiency of the proposed method， an experiment is conducted on the image sequences named as “Nato_camp”， “Tree”， and “Dune” from the TNO Image Fusion Dataset ^［34］. Table 4 shows the average processing time of all methods for a frame. Compared with the DWT， DTCWT， SWT， and WPT methods， the proposed method is more time-consuming because these four methods contain one types of multi-scale decomposition while the proposed method contains two， i.e.， the multi-scale decomposition using multi-scale Gaussian filtering and the multi-scale morphological decomposition， as mentioned in Sec.2. Compared with the NSCT and NSST methods which also contain two kinds of multi-scale decomposition， the proposed method is far more efficient mainly because the design of the multi-directional filter banks for NSCT and NSST is relatively complex and the processing speed of multi-directional filtering is much lower than that of multi-scale morphological operations.

Experiments on both visual quality and objective assessment demonstrate that although adopting the simple AVG-ABS rule， the proposed method does not generate certain artifacts or distortions and performs very well in aspects like information preservation and contrast improvement. Under the premise of ensuring image fusion quality， the proposed method is also proved computationally efficient. The proposed method provides an option for the fusion situations needing both high quality and particularly computational efficiency， such as fast high-resolution images fusion and video fusion.