(8)

CrossMark

Available online at www.sciencedirect.com

ScienceDirect

Procedia Computer Science 57 (2015) 554 - 562

3rd International Conference on Recent Trends in Computing 2015 (ICRTC-2015)

Remote sensing image fusion using multithreshold Otsu method in

shearlet domain

Biswajit Biswasa, Kashi Nath Deyb, Amlan Chakrabartia,b*

a University of Calcutta, 92-A.P.C road, Kolkata and Kol-700009, India b University of Calcutta, 92-A.P.C road, Kolkata and Kol-700009, India c University of Calcutta, 92-A.P.C road, Kolkata and Kol-700009, India

Abstract

In remote sensing image fusion, preservation of spectral information and enhancement of spatial resolution are key issues. In this paper, a novel approach of remote sensing satellite image fusion method have been proposed based on Otsu's Multi-thresholding Method (MOM) in shearlet domain. We make that happened in two folds, i) shearlet transform (ST) is applied in each high-spatial-resolution Panchromatic (PAN) and multi-spectral (MS) image separately, ii) the updated low frequency sub-band shearlet coefficients from decomposed shealet images are composed by the MOM method and select largest low-pass band automatically. The process of different high-pass sub-band shearlet coefficients have been discussed in detail. For obtaining the fused result we use the inverse shearlet transformation (IST). The experimental results show that the proposed method outperforms many state-of-the-art techniques in performance assessment.

© 2015TheAuthors.Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer-review under responsibility of organizing committee of the 3rd International Conference on Recent Trends in Computing 2015 (ICRTC-2015)

Keywords: Remote sensing image fusion; shearlet transform; multithreshold; mutual information;

1. Introduction

Presently the remote sensing satellite image fusion is a promising research field, in the area of the satellite imaging system. Using remote sensing images from multiple satellite image sensors increases robustness and enhances accuracy in remote sensing research1. However, the image fusion can broadly defined as the combination of visual information contained in any number of input images into a single fused image without introducing distortion or information loss2,3. At the pixel-level of the image fusion, the simplest method is to average the input images to produce a fused version. This method is seriously degrade the brightness of the input images4. The pixel-level image fusion include several schemes, such as Intensity-Hue-Saturation (IHS) method1,7, Principal components analysis (PCA) method1,7, Bravery Transform (BT)16, Gram-Schmidt technique (GST)13, Laplacian pyramid method (LP)18, various statistical4,5 and soft computing approach1,2,6 etc.

* Corresponding author. Tel.: +0-000-000-0000 ; fax: +0-000-000-0000. E-mail address: biswajit.cu.08@gmail.com

1877-0509 © 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer-review under responsibility of organizing committee of the 3rd International Conference on Recent Trends in Computing 2015 (ICRTC-2015) doi: 10.1016/j.procs.2015.07.388

Nowadays, the multi-scale decomposition (MSD) based image fusion methods have been widely reported, for example, the fusion of the multi-focus image in1,3, the fusion of the remote sensing image in9,16 and biological image fusion in3,4,6. Recently, a theory of multidimensional data known as multi-scale geometric analysis (MGA) has been developed. The MGA tools were proposed, namely ridgelet, curvelet, contourlet, etc.9,10,11. The new MGA tools provide higher directional sensitivity more than wavelets10,11. Shearlets, a new multi-scale geometric analysis method that equipped with a rich mathematical structure similar to wavelets, that are related to a multi-resolution analysis. The shearlets form a tight Parseval frame at various scales and directions, and are optimally sparse in representing images with edges. Characteristically, the curvelets are the similar with shearlets11,16. The decomposition of shearlets is more like contourlets, as the contourlet transform consists of an application of the Laplacian pyramid, followed by directional filtering, but for shearlets, the directional filtering is replaced by a shear matrix. An important advantage of the shearlet transform over the contourlet transform is that there are no restrictions on the direction numbers10,16.

The shearlet coefficients are highly beneficial in different scales and direction. Here a new remote sensing image fusion algorithm based on MOM in shearlet domain have been proposed. In this paper, a latest mathematical model for commutation and largest shearlet coefficients selection is done. Furthermore, Singular Value Decomposition (SVD) method is a popular method for feature extraction and used for image fusion14. Conventionally, SVD is a mathematical that transform a number of correlated variables into a several orthogonal variables15. The merits of SDV have utilized used in this method as to reduced over dominate intensity in high spatial panchromatic (PAN) image to avoid the spectral deformation. The algorithm evaluates the corresponding high-pass sub-bands of the images at different decomposition levels and also the directional sub-bands of shearlet transform are evaluates. Finally, most fittest shearlet subbands are process based on MOM and reconstruct fused image by inverse shearlet transform. Experimental results demonstrate that the proposed algorithm offers an efficient fusion technique that means to allow more correct analysis of remote sensing satellite image.

The rest of this paper is organized as follows. Brief introduction of ST is given in Section 2, techniques for construction and modification of low-pass and high-pass band is given in Section 3. The fused results and discussions are briefed in Section 4. The final conclusion is summarized in Section 5.

2. Shearlet transform

Conventionally, the shearlet transform (ST) have been developed based on an affine system with composite dila-tions11,16. Let us briefly discussed the continuous and discrete shearlet transform at the fixed resolution level j as follows:

2.1. Continuous shearlet system

In dimension n = 2, the affine systems with composite dilations are the collections of the form16 is given as: VpQ (f) = [f jM (X) (1)

= |det P|j/2f (QlPjX - k) : j, l e Z, k e Z2}

where, P, Q are 2 x 2 invertible matrices, f e L2(R2), mostly represented as follows:

J]l(f,f jk.i>l2 = If"2 (2)

The elements of this system are known as composite wavelets if x¥PQ(f) forms a Parseval frame or tight frame for f e L2(R2). In this system, the dilations matrices Pj is the scale transformations, while the matrices Ql are the area-preserving geometric transformations, such as rotations and shear. As stated in the continuous shearlet transform, we have9,16,

fa,s,k (x) = a-3/4f U-1 V-1 (x - k)) (3)

Where Va = j 0 ^)' U = (Q 0 ), f e L2 (R2), also satisfy the following conditions :

(a) (b) (c)

Fig. 1. The structure of the frequency tiling by the shearlet: (a) The tiling of the frequency plane R2 induced by the shearlet;(b) The size of the frequency support of a shearlet (c) The multi-scale (k = 2) and multi-directional (l = 2) decompositions of ST.

1. f(e) = f(ei, £2) = fi(ei)f2teM);

2. tf/1 e C°°(R) supp c [-2, -1/2] U[1/2,2], where, is continuous wavelet;

3. ^ e C°°(R) supp ^ c [-1,1], ^ > 0 But H^ll = 1

For i^a,s,k, a e R+, s e R+ and k e R+, for any f e L2(R2), is called shearlet 16, more specifically a collection of wavelets with different scales. Here, the anisotropic expansion matrix Us is associate with the scale transform and the shear matrix Va specify the geometric transformation. Mostly, a = 4 and s = 1. Where a, s, k are specify as scale transformations, shear direction and translation vector respectively.

2.2. The discrete shearlet system

The process of the discrete shear transformation can be divided into two steps such as multi-scale subdivision and direction localization9,11,16. Fig. 1 shown the example of decomposition process by the shearlets system.

In this system, at any scale j, let f e L(Z2N). Firstly, the Laplacian pyramid method is used to decompose an image faj-1 into low-pass image fj and a high-pass image fj with Nj = Nj-1/4, where faj-1 e L(1?N, i), fj e L(Z2N) and fj e L(ZN i). After decomposition, estimated f b on a pseudo-polar grid with the one-dimensional band pass filter

based on the signal components, which generate a matrix Dfj. Then, we apply a band-pass filter on matrix Dfb to reconstruct the Cartesian sampled values straightforwardly and also apply the inverse two-dimensional Fast Fourier Transform (FFT) to reconstruct the image16.

2.3. Singular value decomposition ( SVD)

Singular value decomposition (SVD) is arbitrarily decomposed into a singular value matrix with contain only a few non-zero values, this technique can be use to feature extraction from an image14,15. Conventionally, SVD of an m x n matrix A is given by14,15:

A = UxEAVA (4)

where the columns of the m x n matrix UA are known as left singular vectors, the rows of the n x n matrix VA contain the elements of the right singular vectors, and the diagonal elements of the n x n diagonal matrix HA = diag(^1, ••• , iTn) are called as the singular values. Besides, > 0 for 0 < i < q and = 0 if (q + 1) < i < n. By the rule, the arranging of the singular vectors is determined by decreasing order of singular values, with the highest singular value existing in the upper left corner of the HA matrix, specifically, (i1 > i2 > • • • > iq > 0). The concept and method of SVD have used in fusion scheme as follows:

• First, PAN low-pass sub-band image is transformed by SVD transform and achieved the singular values with corresponding upper and lower triangular matrices and singular value component of individual band.

• Second, each singular value is reduced by a weighted factor.

• Finally, PAN low-pass sub-band image obtained by inversion SVD.

2.4. Multilevel Otsu's thresholding technique

In common, the Otsu's multilevel thresholding method (MOM) is most popular threshold selection method for image segmentation and shape measures. However, it is an exhaustive search to estimate the criterion for maximizing the between-class variance12,13. MOM method select a global threshold value by maximizing the differentiability of the cluster in the gray levels12,13.

Let a 2D image I with L gray levels (0,1, ••• , L - 1) and the number of pixel at level l e L is hl in the histogram. Thus, the total number of pixels N be N = h0 + h1 + ••• + hL-1. Generally, from the histogram of a gray level image, the occurrence probability of a pixel in gray level l is given as follows:

Pl = N Pl * 0 X Pl = 1 (5)

In order to segment a 2D image in to m number of clusters (c0, c1 ••• , cm-1 e cm), then m - 1 number thresholds such as (t0, t1 ••• , Tm-2 e rm-1) must be selected. Thus, the cumulative probability pm and average gray level pm for each cluster cm are given by:

Pm = Y Pl and pm = Y l • pM (6)

lecm iecm \Pm'

Hence, the average intensity pT and between-class variance < of image I are estimated by:

L-1 m-1 m-1 m-1

Pr = X l ^ Pl = X PmPm and <2 = ^ Pm • (Pm - Prf = ^ Pm^m - £ (7)

l=0 m=0 l=0 m=0

Thus, the optimal thresholds (r0, ri[ ••• , r"m_2) is determined by maximizing the between-class variance by:

{r0,r ••• ,T*m_2) = argmax {(t0,T1 ••• ,Tm-2)} (8)

0<T0<...<Tm-1<L

The Multilevel Otsu's thresholding technique is used in the proposed fusion scheme as follows:

• The MOM is applied on the both PAN and MS images and obtained the corresponding multiple threshold within a vector based on a given clusters number m. In this paper, we considered m = 5.

• From each multiple threshold vector, we are estimate the average thresholds and denoted by a and y separately.

• Finally, a and y are used to define an adaptive factor to modify the low-pass sub-bands of MS and PAN.

3. The proposed fusion algorithm

In this paper, the proposed fusion algorithm named as MOMST. The main steps of the MOMST is shown in Fig. 2. The overall process is introduced as follows:

3.1. ComPOsition of low-Pass sub-bands

Low-pass sub-bands of the ST give the approximation of original images9,16. However, to preserve spectral information and excellent image composition, we are evaluate shearlet low-pass sub-bands by MOM and finally select largest low-pass sub-bands successfully.

Let Ls(i, j) is the low-pass sub-band located at (i, j), s = PAN, MS. In this paper, LPAN consider as high spatial image (PAN). Thus, before integration the LPAN is reduced by SVD as follows:

[ Upan, Zpan, VPAN] = SVD (Lpan) (9)

where each YPAN singular value is divided by a weighted factor W, that denoted as £PAN = ^W * log(W)3/2(£PAN./ VW), where W = r * c, r and c are row and columns of PAN. Finally, we obtained the reduce LPAN as LPAN = VPAN * £PAN *

Fused image

Fig. 2. Main steps of MOMST fusion algorithm. However, the fusion of low-pass sub-bands LF(x, y), which are represented as follows:

lf(U j) = {Lms (i, j) * n if Lms (i, j) > Lpan(i, j), otherwise Lpan(i, j) * n (10)

where n = V1 + a2 + y2 and a, y are parameters that are estimated by MOM method.

Before MOM method, we use gradient operator on both LPAN, LMS matrices to estimate the over all horizontal and vertical variation. The gradient operator of an image with size N x N is given by1:

x = (AL)x, = [^ - L"j,if'' < N , y = (AL? . = L1 j - LiJ,if< N (11)

V h,] [0, if i = N y K >l'] (0, if i = N K '

where x and y are horizontal and vertical image gradient of the image respectively.

To evaluate the two parameters a and y, we have successfully used the advantage MOM. However, the core computation can be summed up as follows:

• Foremost, we compute gradient matrices Xh, Yh, Xv, and Yh, from low-pass sub-band matrix LX and LY by using Eqs.11 and11. Where h, v stand for horizontal and vertical direction of LX and LY respectively.

• Next, we estimate absolute disparity from each horizontal and vertical matrix Xh, Yh, Xv, and Yh, and expressed as Dp = (|Xh - Yh|) and Dvp = (|Xv - Yv|) separately.

• After disparity measure, MOM is used on Dhp, Dvp. We estimated two parameters a and y from Dhp and Dvp and these parameter are used to define the factor n. So the a and y are estimated by:

/log(ZGupAN)) ZpAN , „ (logZ(MMS)) ZmS

£pan = V-, a -rr~ and ZMS = V-, Y -—— (12)

v m Zpan + Zms y m Zpan + Zms

Where ¡iPAN, nMS are the multilevel threshold vectors of PAN and MS and m is number of clusters.

3.2. Composition of high pass sub-bands

High-pass sub-bands of the ST provide the details information of the image, such as edge corner, etc.9,11,16. In MOMST, a new decision map have been developed to upgrade the high-pass sub-bands and an efficient fusion rule that combined new high-pass sub-bands automatically.

Let Hjk(x, y) be the high-pass coefficient at the location (x,y) in the lth sub-band at the kth level, s = PAN, MS. In short, the process of high-pass sub-bands of MOMST is discussed as follows:

For current sub-band Hlsk, let Qs,h is the total of H1^ and the other horizontal sub-bands Hm,k in same level k. It can be estimated by:

Qsp = Zi>s-1,k - rtfy (13)

k=\ /=1

Similarly, let Qsv is the sum of Hl.;k and all the other vertical high-pass sub-bands Hm'n in the different finer levels. It can be evaluated by:

Qsv =ZZl(HS,k - H7'n)\ (14)

Cs, v '

l,m=1 k,n=1

To determine the parameter ps,h along horizontal and the parameter ps,v along vertical high-pass sub-bands for the

present high-pass sub-band H,k, we perform:

_ Qs,h _ Qs,v rx

ps,h = (Qsh + Qsv) ps,v = (Qs,h + Qsv) ()

The parameter ps,h is a relationship between Hlsk and other neighbor high-pass sub-bands in the same horizontal plane, and ps,v relationship between H'sk in the different vertical plane with its corresponding neighbor sub-bands. Finally, to compute the new coefficients H'!knew using by parameter ps,h and ps v, we compute:

Hlsknew = Hlsk X ^ 1 + plh + p?,v (16)

Finally, the fused shearlet coefficients Hlp(x, y) is obtained by the following computation:

HF(x, y) = [h1* (i, j), if HlPk > Hl* otherwise Hl* (i, j) (17)

* I PAN,newK ' PAN,new MS,new MS,newv ' v '

The overall procedure of the MOMST fusion algorithm can be summarized as follows:

• The PAN and MS images must be registered.

• These source images are decomposed into shearlet images by ST.

• Best low-pass sub-bands are estimated by MOM and fused by using Eq. 10.

• Finest high-pass shearlet coefficients are modified and select largest shearlet coefficients by using Eqs. 16, 17.

• Largest directional shearlet coefficients are selected.

• The fused image is reconstructed by 1ST.

4. Experimental Results and Discussion

In this section, we conducted two illustrative examples to demonstrate the efficiency of our proposed method MOMST. Firstly, the image select as QuickBird satellite image with size (400 x 400). The MS image of QuickBird data has four bands: blue band (450 - 520 nm), green band (520-600 nm), red band (630 - 690 nm) and near infrared (NIR) band (760 - 900 nm), the resolution of MS is 2.44 m. The resolution of PAN image (450 - 900 nm) is 0.70 m19. The simulations are performed in MATLAB R2013a on PC with Intel(R) Core(TM)-2 Quad CPU/2.33G/6G. To evaluate the performance of MOMST fusion in PAN and MS image, a valid experiment is implemented and the results are compared in term of qualitatively and quantitatively. In order to validate of the experiment, we compared the performance of the proposed MOMST method with six different fusion schemes, i.e. (1) Intensity-Hue-Saturation (IHS)1,7?, (2) Principal component analysis (PCA)1,7 (3) Laplacian pyramid Technique (LPT)1,8 (4) Discrete wavelet transform (DWT)3, (5) Bravery Transform (BT)1A6, (6) Gram-Schmidt technique (GST) 1,3,5respectively.

In performance assessment, four popular metrics such as Entropy (EN)1,n, Average Gradient (AG)1,n, Mutual Information (MI)17 and QAB/F 17,18 are considers. The results by EN, AG, MI and QAB/F are shown in Table 1 and Table 2 of Figs. 3 and 4 separately. Note that the highest score in each row of Table 2 are shown in bold. From that table, it can be noted that the MOMST method outperforms than other methods in term of evaluation metrics. Moreover, we also perform a visual comparison on image set set-I as shown in Fig. 3, set-II as shown in Fig. 4. In Figs. 3-4 are demonstrated the resultant fused images obtained from the IHS, PCA, LP, DWT, BT, GST and MOMST respectively. Experiments illustrated that MOMST approach can resolve spectral distortion problems and successfully preserve the spatial information of fuse image. Therefore, the fused image obtained from the proposed method MOMST provide higher quality fusion image than the other methods in state-of-the-art techniques. It can be concluded that the MOMST is more efficient and stable.

. | . л - p E _______t.

(d) PAN (e) MS (f) MOMST

Fig. 3. (a) Multispectral QuickBird image (QuickBird, 4-m resolution, 400 x 400 pixels);(b) Panchromatic QuickBird image (QuickBird, 1-m resolution, 400 x 400 pixels);(c) MOMST.

Table 1. Statistical result analysis of MOMST.

Image qF/ab EN MI AG

Data I 0.3413 5.1271 4.2138 40.3741

Data II 0.4139 5.2198 4.1547 40.9529

Table 2. Quantitative comparison and result analysis.

Method BT GST IHS PCA LPT DWT MOMST

QF/AB 0.3306 0.3468 0.3526 0.3548 0.3597 0.3694 0.3863

EN 4.4621 4.4969 4.5219 4.5837 4.6970 4.7168 4.9308

MI 3.4614 3.4720 3.5097 3.5176 3.5921 3.6091 3.9376

AG 37.372 37.259 37.589 37.764 37.984 38.619 41.262

5. Conclusion

In this paper we have visualized and successfully implemented, the remote sensing satellite image fusion method based Otsu's multilevel thresholding method on shearlet domain. Most importantly we are able to preserve color information as well as maintaining the spectral resolution of the fused image. As our proposed MOMST technique for image fusion compered with other available state-of-art techniques, where MOMST tenders better result.

(g) DWT (h) LPT (i) MOMST

Fig. 4. (a) Multispectral QuickBird image (QuickBird, 4-m resolution, 400 X400 pixels after resampling and resizing);(b) Panchromatic QuickBird image (QuickBird, 1-m resolution, 400 X 400 pixels after resizing) image fusion by various techniques;(c) GST; (d) IHS;(e) PCA;(f) BT; (g) DWT; (h) LPT; (i) MOMST.

References

1. Konstantinos N G, Comparison of nine fusion techniques for very high resolution data, Photogramm. Eng. Remote Sens. (2008) 74(5) pp.647659.

2. Zhijun W, Ziou D, Armenakis C, Deren Li, Qingquan Li, A comparative analysis of image fusion methods, Geoscience and Remote Sensing, IEEE Transactions on, 2005 43(6) pp.1391-1402.

3. Pajares G, Cruz JM, A wavelet-based image fusion tutorial, Pattern Recognition, (2004) 37(9) pp.1855-1872.

4. Yang B, Jing Z, Zhao H, Review of pixel-level image fusion, Journal of Shanghai Jiaotong University (Science), (2010) 15(1) pp.6-12.

5. Pohl C, Van G J L, Review article Multisensor image fusion in remote sensing: Concepts, methods and applications, IJRS,( 1998) 19(5) pp.823-854.

6. Simone G, Farina A, Morabito FC, Serpico SB, Bruzzone L, Image fusion techniques for remote sensing applications, INFORMATION FUSION, (2002) 3 pp.3-15.

7. He C, Liu Q, Li H, Wang H, Multimodal medical image fusion based on IHS and PCA, Procedia Engineering, (2010) 17 pp.280-285.

8. Thevenaz P, Unser M, A pyramid approach to sub-pixel image fusion based on mutual information, (1996) (1) pp.265-268.

9. Wang H, Liu Y, Xu S, An image fusion algorithm based on shearlet, (2013) pp.1625-1628.

10. Do M N, Vetterli M, The contourlet transform: an efficient directional multiresolution image representation, IEEE Trans. Image Process. (2005) 14 pp.2091-2106.

11. Deng C, Wang S, ChenX, Remote Sensing Images Fusion Algorithm Based on Shearlet Transform, (2009) 3 pp.451-454.

12. Liao PS, Chen TS, Chung PC, A Fast Algorithm for Multilevel Thresholding, JISE, (2001) 17 pp.713-727.

13. Huang DY, Lin TW, Hu WC, A Fast Algorithm for Multilevel Thresholding, International journal of innovative computing, information and control, (2011) 7(10) pp.5631-5644.

14. Pinkus A R, etal, Investigation of image fusion procedures using optimal registration and SVD algorithms, NAECON (2009) pp. 231-235.

15. Nasir H, etal, Singular value decomposition based fusion for super-resolution image reconstruction,Signal Processing-Image Communication, (2011).

16. Easley G, Labate D, Lim WQ, Sparse Directional Image Representations using the Discrete Shearlet Transform, Applied Computational Harmonic Analysis, (2008) 25 pp.25-46.

17. Xydeas C, Petrovic V, Objective image fusion performance measure, Electronic Letters. (2000), 36(4) pp.308-309.

18. Chen, Yin, Xue, Zhiyun, Blum, Rick S, Theoretical analysis of an information-based quality measure for image fusion, Inf. Fusion, (2008) 9(2) pp.161-175.

19. http://glcf.umiacs.umd.edu/data.