- Email: [email protected]
Discrete multivariate gray model based boundary extension for bi-dimensional empirical mode decomposition
Signal Processing 93 (2013) 124–138
Contents lists available at SciVerse ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Discrete multivariate gray model based boundary extension for bi-dimensional empirical mode decomposition Zhi He, Qiang Wang, Yi Shen n, Yan Wang School of Astronautics, Harbin Institute of Technology (HIT), Harbin 150001, China
a r t i c l e in f o
abstract
Article history: Received 9 December 2011 Received in revised form 24 May 2012 Accepted 9 July 2012 Available online 20 July 2012
The bi-dimensional empirical mode decomposition (BEMD) has attracted extensive attention recently by virtue of its high performance in adaptive image processing. Unfortunately, this promising technique does not necessarily yield fruitful results due to the boundary effects. Motivated by the discrete multivariate gray model, we propose a boundary extension framework for mitigating the boundary effects of BEMD. In greater detail, followed by verifying the equivalence between the continuous and discrete multivariate gray model theoretically, a first-order three-variable discrete multivariate gray model D-GMC(1,3), which is derived from the continuous multivariate gray model with convolution integral C-GMC(1,N), is utilized to predict the middle pixels of each extended block in terms of existing border. Specifically, the coordinates and pixels of the image are respectively regarded as relative data series and characteristic data series of D-GMC(1,3). Experimental results on a set of widely used images indicate that the proposed approach can achieve qualitative and quantitative improvements within appropriate processing time by comparing with other three generally acknowledged methods, i.e. the original BEMD, symmetrical extension as well as texture synthesis based BEMD. & 2012 Elsevier B.V. All rights reserved.
Keywords: Bi-dimensional empirical mode decomposition (BEMD) Boundary effects Discrete multivariate gray model
1. Introduction Over the last decades, the bi-dimensional empirical mode decomposition (BEMD) has become one of the most significant image processing techniques owing to its adaptive nature and versatile capability. In analogy to the one-dimensional (1-D) empirical mode decomposition (EMD) [1–3], the BEMD is also adaptive, empirical and intuitive due to its fully data-driven property. To date, the development of BEMD is twofold: theory and application. More precisely, Nunes et al. developed the EMD into two-dimensional (2-D) version, namely BEMD [4], in which the morphological reconstruction and radial basis
n Corresponding author. Tel.: þ 86 451 86413411; fax: þ86 451 86418378. E-mail address: [email protected] (Y. Shen).
0165-1684/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sigpro.2012.07.009
function (RBF) are implemented to detect the image extrema and calculate the surface interpolation, respectively. Bhuiyan et al. evaluated the existing surface interpolation techniques and concluded that the RBF based interpolators show promising performance [5]. To alleviate the scale mixing, a noise-assisted method termed as multi-dimensional ensemble EMD (MEEMD) was introduced by Wu et al. [6]. Recently, Yeh established a complex extension of BEMD, which took advantage of the perfect reconstruction of the four quadrant spectra and suitable for 2-D complex signal processing [7]. Fleureau et al. studied the multivariate EMD in terms of redefining the local mean operator [8]. Moreover, a number of fast algorithms based on finite elements [9], order-statistics filter [10] or adaptive spline-wavelets [11] have also been exploited. On the other hand, different kinds of applications are booming up nowadays. For instance, Linderhed [12,13] and Tian et al. [14] developed
Z. He et al. / Signal Processing 93 (2013) 124–138
a totally different image compression approach inspired by the multi-resolution properties of bi-dimensional intrinsic mode function (BIMF). A BEMD-morphology approach (BMA), which was proposed by Chen et al., is effective in handling non-uniform illumination in bridge coating assessment [15]. Huang et al. applied the BEMD to decompose the gravity data for the Tongshi gold field and obtained a clearer spatial distribution than that of Fourier transform [16]. As stated in [17], the desired edges of image were detected by a fast and adaptive BEMD algorithm abbreviated as FABEMD. Though BEMD has recently developed a valuable tool for image processing, how to mitigate boundary effects (see Fig. 1) remains an open problem, generated by the lack of extrema on the boundary or beyond the original range of the image. Indeed, the extrema may be sparse thus the envelopes detected by interpolating local extrema cannot cover the existing region of the image. Worse still, those drawbacks may propagate into the interior part of the image and ‘‘pollute’’ the full decomposition results. Much work has been carried out in the literature to handle this issue, whose merit is to extend data explicitly or implicitly out of the original region as suggested by Huang et al. in 1-D case [1]. Typically, the
125
symmetrical or periodic extension based BEMD extends the given image symmetrically or periodically. Nevertheless, the performance of this method is far from satisfactory in case the image is asymmetrical or nonperiodic. Moreover, Liu et al. introduced an extension technique inspired by nonparametric sampling based texture synthesis method [18,19]. Unfortunately, two main problems related to this approach: (1) it requires that the given image is extended accurately enough before decomposition; (2) it desires quantitative judgment of various forms of boundary processing. Other than the above-mentioned methods, we propose an alternative boundary extension method motivated by the discrete multivariate gray model [20–22]. Indeed, as a valuable forecasting technique, the gray model has successfully alleviated the end effects of Hilbert–Huang transform in our previous work [23]. More critically, in this paper, followed by evaluating the equivalence between the continuous and discrete multivariate gray model theoretically, a first-order three-variable discrete multivariate gray model termed as D-GMC(1,3), which is derived from the continuous multivariate gray model with convolution integral C-GMC(1,N), is associated with the extension of the given image. In general, we regard the coordinates and pixels of the image as relative data series and characteristic data series of D-GMC(1,3), respectively. Once the boundary of the given image is extended by D-GMC(1,3), BEMD is devoted to extracting the essential components, namely BIMFs. It is noteworthy that the bi-dimensional index of orthogonality (BIO) is introduced to evaluate various boundary processing methods quantitatively. The outline of the paper is organized as follows: Section 2 introduces the principle of the multivariate gray model. This is followed by a detailed description of our proposed boundary extension scheme in Section 3. In Section 4, comparative analysis of the widespread methods is carried out in terms of extensive experiments on standard images. The conclusion is drawn in Section 5. 2. The multivariate gray model 2.1. The continuous multivariate gray model
Fig. 1. An example of boundary effects during extracting the first BIMF. (a) The 3-D graphic of ‘‘LENA’’ depicted in Fig. 6(a). (b) The mean envelope in the process of obtaining the first BIMF.
As mentioned earlier, the boundary effects of BEMD are mainly due to lacking of extrema on the boundary or beyond the original region of the image. To alleviate this annoying deficiency, a boundary extension algorithm inspired by the discrete multivariate gray mode is proposed in this paper. Alternatively, the keynote of this section focuses on the overview of continuous multivariate gray model with convolution integral, termed as C-GMC(1,N) [21], which is recently proposed by Tien in view of the incorrect solution of the original gray model GM(1,N). Indeed, it is the foundation of discrete multivariate gray model as will be described in the following section. Assume the original characteristic data series and ð0Þ ð0Þ ð0Þ relative data series are Xð0Þ 1 ¼ fx1 ð1Þ,x1 ð2Þ, . . . ,x1 ðnÞg ð0Þ ð0Þ ð0Þ ð0Þ and Xi ¼ fxi ð1Þ,xi ð2Þ, . . . ,xi ðn þmÞg, i ¼ 2,3, . . . ,N, respectively, where m is the number of entries to
126
Z. He et al. / Signal Processing 93 (2013) 124–138
be predicted, xð0Þ and xð0Þ ðkÞ Z0, 1 ðkÞ Z 0,k ¼ 1,2, . . . ,n i ð0Þ k ¼ 1,2, . . . ,n þ m, i ¼ 2,3, . . . ,N. xi ðkÞ can be considered as integer values sampled from a continuous function. Pk ð0Þ ð1Þ Given xð1Þ l ¼ 1 x1 ðlÞ, k ¼ 1,2, . . . ,n and xi ðkÞ ¼ 1 ðkÞ ¼ Pk ð0Þ k ¼ 1,2, . . . ,n þ m, i ¼ 2,3, . . . ,N as the l ¼ 1 xi ðlÞ, first-order accumulated generating operation (1-AGO) of xið0Þ ðkÞ, the whitenization form of C-GMC(1,N) can then be defined as ð1Þ N 1 X dx1 ðtÞ bi xð1Þ ðtÞ þbN , þ axð1Þ 1 ðtÞ ¼ iþ1 dt i¼1
2.2. The discrete multivariate gray model Versatile as C-GMC(1,N), the improvement of prediction accuracy remains a challenging problem. So far, much work has been carried out in the literature to handle this issue [21,22]. Despite the partial improvements, all these methods do not solve the problem intrinsically. As stated in [20], the coefficients are evaluated by discrete equation, whereas the modeling and predicting are given by continuous equations. Consequently, the paradox between discrete estimation and continuous prediction is just the cause of forecast error in C-GMC(1,N). To counter this effect, we introduce the discrete multivariate gray model, in which the parameter estimation, modeling construction together with data predication are achieved by discrete equations and there is no contradiction between discrete estimation and continuous prediction, and therefore, has high prediction accuracy.
t ¼ 1,2, . . . ,n þ m ð1Þ
where a is gray developmental coefficient, bi ,i ¼ 1,2, . . . ,N are associated coefficients and rp ¼0 in Eq. (33) of [21] taking into account the actual situation of this paper. Letting k1 ¼ ½a b1 b2 . . . bN T , we can determine k1 as follows
k1 ¼ ½a b1 b2 . . . bN T ¼ ðBT1 B1 Þ1 BT1 Y 1
Definition 1. Assume the original characteristic data ð0Þ series Xð0Þ 1 and relative data series Xi , i ¼ 2,3, . . . ,N with n entries are the same as those of Section 2.1. The
ð2Þ
ð0Þ T where Y 1 ¼ ½x1ð0Þ ð2Þ xð0Þ 1 ð3Þ . . . x1 ðnÞ and
2
ð1Þ 12 ðxð1Þ 1 ð1Þ þ x1 ð2ÞÞ 6 6 1 ðxð1Þ ð2Þ þ xð1Þ ð3ÞÞ 6 1 2 1 B1 ¼ 6 6 ^ 4 ð1Þ 12 ðxð1Þ 1 ðn1Þ þx1 ðnÞÞ
ð1Þ 1 ð1Þ 2 ðx2 ð1Þ þ x2 ð2ÞÞ ð1Þ 1 ð1Þ 2 ðx2 ð2Þ þ x2 ð3ÞÞ
^
&
^
ð1Þ 1 ð1Þ 2 ðx2 ðn1Þ þx2 ðnÞÞ
ð1Þ 1 ð1Þ 2 ðxN ðn1Þ þxN ðnÞÞ
aðt1Þ ¼ xð1Þ þ 1 ð1Þe
aðt1Þ ¼ xð1Þ þ 1 ð1Þe
Z
Z
t
eaðttÞ
1 t
N 1 X
! bi xð1Þ ðtÞ þ bN iþ1
ð1Þ xð1Þ 1 ðkÞ ¼ a1 x1 ðk1Þ þ
dt
Z
1
t
eaðttÞ
1
bN ð1eaðt1Þ Þþ a
Z
t
eaðttÞ
1
N 1 X i¼1
N 1 X
bi xð1Þ ðtÞ dt iþ1
i¼1
ð3Þ
þ 1
aðktÞ
e
bN ð1eaðk1Þ Þ a
N 1 X
bi xð1Þ i þ 1ð
tÞ dt
ð4Þ
i¼1
ð1Þ
where x^ 1 ðkÞ denotes the predicted value of xð1Þ 1 ðkÞ. Therefore, the predicted value of x1ð0Þ ðkÞ, k ¼ 2,3, . . . ,n þ m yields ð0Þ ð1Þ ð1Þ x^ 1 ðkÞ ¼ x^ 1 ðkÞx^ 1 ðk1Þ
bi 2
ðxð1Þ ðk1Þ þ xð1Þ ðkÞÞ þ bN iþ1 iþ1 ð6Þ
bi xð1Þ ðtÞ dt iþ1
P ð1Þ where 1 eaðttÞ ð N1 i ¼ 1 bi xi þ 1 ðtÞ þ bN Þ dt is just the convolution integral that the C-GMC(1,N) differs from the original GM(1,N). Consequently, the predicted series of xð1Þ 1 ðkÞ,k ¼ 2,3, . . . ,n þ m are deduced from
k
N 1 X i¼1
Rt
Z
3
7 17 7 7 ^7 5 1
i¼1
eaðttÞ bN dt þ
ð1Þ aðk1Þ þ x^ 1 ðkÞ ¼ xð1Þ 1 ð1Þe
1
first-order N-variable discrete multivariate gray model, namely D-GMC(1,N), can be defined as
Solving Eq. (1) leads to ð1Þ aðt1Þ xð1Þ þ 1 ðtÞ ¼ x1 ð1Þe
ð1Þ 1 ð1Þ 2 ðxN ð1Þ þ xN ð2ÞÞ ð1Þ 1 ð1Þ 2 ðxN ð2Þ þ xN ð3ÞÞ
ð5Þ
Pk
ð0Þ where xð1Þ l ¼ 1 xi ðlÞ, k ¼ 1,2, . . . ,n þ m, i ¼ 1,2, . . . ,N i ðkÞ ¼ and k2 ¼ ½a1 b1 b2 . . . bN T is the coefficient.
In analogy to C-GMC(1,N), the coefficient k2 is evaluated as
k2 ¼ ½a1 b1 b2 . . . bN T ¼ ðBT2 B2 Þ1 BT2 Y 2 where
Y 2 ¼ ½xð1Þ 1 ð2Þ
2
xð1Þ 1 ð1Þ 6 6 xð1Þ ð2Þ 6 1 B2 ¼ 6 6 ^ 4 ð1Þ x1 ðn1Þ
xð1Þ 1 ð3Þ
...
T xð1Þ 1 ðnÞ
ð1Þ 1 ð1Þ 2 ðx2 ð1Þ þ x2 ð2ÞÞ ð1Þ 1 ð1Þ ðx ð2Þ þ x 2 ð3ÞÞ 2 2
^ ð1Þ 1 ð1Þ 2 ðx2 ðn1Þþ x2 ðnÞÞ
and ð1Þ 1 ð1Þ 2 ðxN ð1Þ þ xN ð2ÞÞ ð1Þ 1 ð1Þ ðx ð2Þ þ x N ð3ÞÞ 2 N
&
ð7Þ
^ ð1Þ 1 ð1Þ 2 ðxN ðn1Þ þ xN ðnÞÞ
1
3
7 17 7 7 ^7 5
1
Given the numerous coefficients of D-GMC(1,N), the following intermediate matrix U is introduced to avoid excessive matrix computations: 2
u11 6 u 6 21 U¼6 6 ^ 4 uN þ 1,1
u12
u22
^
&
uN þ 1,2
u1,N þ 1
3
7 u2,N þ 1 7 7 ¼ ðBT B2 Þ1 2 7 ^ 5
uN þ 1,N þ 1
ð8Þ
Z. He et al. / Signal Processing 93 (2013) 124–138
Substituting Eq. (8) into Eq. (7) leads to 2 3 2 3 2 a1 xð1Þ u12 u1,N þ 1 u11 1 ð1Þ 6b 7 6 17 6 7 6 6 u22 u2,N þ 1 7 6 xð1Þ 6 7 6 u21 1 ð2Þ 76 b 7 6 k2 ¼ 6 6 2 7¼6 ^ 7 6 ^ & ^ ^ 6 ^ 7 4 5 4 4 5 uN þ 1,1 uN þ 1,2 uN þ 1,N þ 1 xð1Þ ðn1Þ
bN
1
2
127
ð1Þ 1 ð1Þ 2 ðx2 ð1Þ þ x2 ð2ÞÞ ð1Þ 1 ð1Þ 2 ðx2 ð2Þ þ x2 ð3ÞÞ
^
&
^
ð1Þ 1 ð1Þ 2 ðx2 ðn1Þ þ x2 ðnÞÞ
ð1Þ 1 ð1Þ 2 ðxN ðn1Þ þxN ðnÞÞ
ð1Þ 1 ð1Þ 2 ðxN ð1Þ þ xN ð2ÞÞ ð1Þ 1 ð1Þ 2 ðxN ð2Þ þ xN ð3ÞÞ
3 n X ðxð1Þ ðk1Þ þxð1Þ ðkÞÞxð1Þ ð1Þ 1 ðkÞ i i þu1,N þ 1 u11 xð1Þ u1i xð1Þ 6 7 1 ðk1Þx1 ðkÞ þ 1 ðkÞ 2 7 2 ð1Þ 3 6 i¼2 k¼2 k¼2 k¼2 6 7 x1 ð2Þ 6 7 ð1Þ ð1Þ ð1Þ n N n n X X ð1Þ X X 6 ð1Þ 7 6 7 ðx ðk1Þ þx ðkÞÞx ðkÞ 1 i i 6 x ð3Þ 7 6 7 þu2,N þ 1 x1 ðk1Þxð1Þ u2i xð1Þ u21 6 1 7 6 7 1 ðkÞ þ 1 ðkÞ 2 6 7¼6 7 i¼2 k¼2 k¼2 k¼2 6 ^ 7 6 7 4 5 6 7 ^ 6 7 ð1Þ x1 ðnÞ 6 7 ð1Þ ð1Þ ð1Þ n N n n X X X X 6 7 ðxi ðk1Þ þxi ðkÞÞx1 ðkÞ ð1Þ ð1Þ ð1Þ 4 uN þ 1,1 þuN þ 1,N þ 1 x1 ðk1Þx1 ðkÞ þ uN þ 1,i x1 ðkÞ 5 2 i¼2 k¼2 k¼2 k¼2 n X
N X
Setting 8 n X > ð1Þ > > xð1Þ > CðjÞ ¼ uj1 1 ðk1Þx1 ðkÞ > > > k¼2 > > > > N n < X X ðxð1Þ ðk1Þ þxð1Þ ðkÞÞxð1Þ 1 ðkÞ i i uji DðjÞ ¼ 2 > > i ¼ 2 k ¼ 2 > > > n > X > > > EðjÞ ¼ u xð1Þ > j,N þ 1 1 ðkÞ > : k¼2 the coefficient k2 in Eq. (9) can be rewritten as 2 3 2 3 a1 Cð1Þ þ Dð1Þ þEð1Þ 6b 7 6 17 6 7 Cð2Þ þ Dð2Þ þEð2Þ 6 7 6 7 7¼6 7 b k2 ¼ 6 2 6 7 6 7 ^ 6 ^ 7 4 5 4 5 CðN þ 1Þ þ DðN þ 1Þ þEðN þ 1Þ
n X
ð9Þ
ð10Þ Algorithm 1. Boundary extension algorithm via D-GMC(1,3). 2
ð1Þ x^ ðkÞ ¼ a1 xð1Þ 1 ðk1Þ þ
N 1 X i¼1
bi 2
ðxð1Þ ðk1Þ þ xð1Þ ðkÞÞ þ bN iþ1 iþ1 ð11Þ
and ð0Þ ð1Þ ð1Þ x^ 1 ðkÞ ¼ x^ 1 ðkÞx^ 1 ðk1Þ
ð12Þ
ð0Þ and x^ 1 ðkÞ, k ¼ 2,3, . . . ,n þ m are ð1Þ of x1 ðkÞ and x1ð0Þ ðkÞ, respectively.
where the predicted values We can assess the forecast performance of C-GMC(1,N) or D-GMC(1,N) via the following average relative error: n 1X L nk¼1 k
7 17 7 7 ^7 5 1
Based on the study of the discrete multivariate gray model, we can summarize the proposed boundary extension algorithm in Algorithm 1. In greater detail, 2 the original image IM original ðx,yÞ 2 l ðR2 Þ, ðx,y 2 N þ Þ, whose size is P P, can be predicted by the D-GMC (1,3), in which the coordinates and pixels of the image are respectively regarded as relative data series and characteristic data series. As such, the size of extended image IM extended ðx,yÞ becomes ðP þ 2Q Þ ðP þ2Q Þ with the central P P block being the original image IM original ðx,yÞ.
Having decided the coefficient k2 , the predicted values ð0Þ of xð1Þ 1 ðkÞ and x1 ðkÞ,k ¼ 2,3, . . . ,n þ m can be determined by
L¼
3T
3. The proposed boundary extension algorithm
bN
ð1Þ x^ 1 ðkÞ
1
ð13Þ
where Lk ¼ 9eðkÞ=x1ð0Þ ðkÞ9, and eðkÞ ¼ x1ð0Þ ðkÞx^ ð0Þ 1 ðkÞ, k ¼ 1,2, . . . ,n. Finally, Theorem 1 verifies the the equivalence between C-GMC(1,N) and D-GMC(1,N). Theorem 1. Assume the gray developmental coefficient a of Eq. (1) is small, then C-GMC(1,N) is equivalent to D-GMC(1,N). Proof. The proof of Theorem 1 is given in Appendix A.
Require: The original image IM original ðx,yÞ 2 l ðR2 Þ,ðx,y 2 N þ Þ with size P P; error e1 and e2 ; size of the extended image ðP þ 2Q Þ ðP þ 2Q Þ; size of searching block T T Ensure The extended image IM extended ðx,yÞ 1: %initialization 2: IMprocess ðx,yÞ ¼ IMoriginal ðx,yÞ 3: %extension procedure 4: while IMprocess ðx,yÞ is not fully filled into ðP þ 2Q Þ ðP þ 2Q Þ do 5: Browse through the original image IMoriginal ðx,yÞ in raster scan order 6: Search IM original ðx,yÞ for a set of blocks (sizes are T T) B ¼ fB1 ,B2 , . . . ,Bl g that fulfill the overlap constraints from above and left (only the about overlap is considered in case filling the first column of IM process ðx,yÞ) within the error tolerance e1 7: Pick one block Brandom randomly from B ¼ fB1 ,B2 , . . . ,Bl g, B0 ¼ Brandom 8: Predict the inner pixel values of B0 via D-GMC(1,3) in terms of the surrounding pixel values of Brandom 9: Get the prediction error L of D-GMC(1,3) 10: 11: 12: 13: 14: 15: 16: 17:
if L 4 e2 then B0 ¼ Brandom %in case the prediction accuracy of D-GMC(1,3) is unsatisfactory end if Calculate the error surface S between the newly generated block B0 and the old blocks on the overlap region Identify the minimum cost path along S and set it as the boundary of B0 [24] Tile the modified block B0 onto the corresponding region of IMprocess ðx,yÞ end while IMextended ðx,yÞ ¼ IM process ðx,yÞ
128
Z. He et al. / Signal Processing 93 (2013) 124–138
characteristic data series X ð0Þ 1 and relative data series ð0Þ X ð0Þ and X of D-GMC(1,3) are given by pixel values, 2 3 horizontal as well as vertical coordinates, respectively, of the chosen block B0 . As long as the prediction error is smaller than e2 , these predicted inner pixels will be considered as a partial extended boundary of IM extended ðx,yÞ. 4. It is highlighted from Fig. 3 that the boundary of B0 is evaluated by the minimum cost path [24] along the error surface S. In this way, we can improve the smoothness and cohesive property of overlapping between two small blocks. In greater detail, Fig. 3 outlined all the three potential overlapping scenarios, horizontal, vertical and corner overlaps, respectively.
The following four remarks further clarify Algorithm 1: 1. As depicted in Fig. 2(a), we divide the boundary to be extended into eight small blocks. It is worth underlining that different combinations of those blocks may yield different extension orders, which are typically shown in Fig. 2(b)–(g). Although various combinations of the blocks may influence the boundary extension results, these influences are quite small in general. From this viewpoint, we choose the extension order according to Fig. 2(b) without loss of generality. 2. It is illustrated in Fig. 3 that the extension order of original image is A-B-C-D, which is combined from the eight small blocks of Fig. 2(a). As for ‘‘A’’, it is sequently filled by a few small blocks with sizes T T (see Fig. 3), each of them are randomly picked from a set of candidate blocks. Particularly, the other boundaries, i.e., ‘‘B’’, ‘‘C’’ and ‘‘D’’, can be converted into the same situation as ‘‘A’’ by virtue of folding, rotation as well as transposition. 3. With regard to each small block B0 in Fig. 3, its inner pixels can be predicted by D-GMC(1,3) making full use of all their surrounding pixels. More specifically, the
1
4
6
Having finished illustrating the proposed boundary extension approach, the modified BEMD algorithm inspired by D-GMC(1,3), namely D-BEMD, can be drawn in Fig. 4. It is noteworthy that the main contribution of D-BEMD is twofold:
The original image is extended by D-GMC(1,3) before applying the traditional BEMD.
2
IMoriginal
7
3
5
8
IMoriginal
IMoriginal
IMoriginal
IMoriginal
IMoriginal
IMoriginal
Fig. 2. Schematic division of the extended blocks in Algorithm 1. (a) Eight small blocks to be extended. (b)–(g) Six different combinations of the small blocks in (a).
Z. He et al. / Signal Processing 93 (2013) 124–138
129
Pick a block randomly
B1 B2 Bl
B0
IMoriginal
Enlarged
DMGM (1,3)
D Minimum cost path
A IMoriginal C
B Enlarged
Potential overlapping scenarios
A Horizontal Vertical
Corner
Fig. 3. Schematic presentation of Algorithm 1.
The corresponding parts of the extended BIMFs and final residue obtained from the extended image are extracted as decomposition results of the original image.
4. Experimental results The keynote of this section focuses on assessing the performance of the proposed D-BEMD algorithm by comparing with other three generally acknowledged methods, i.e. the original BEMD, symmetrical extension as well as texture synthesis based BEMD. For the sake of simplicity, the abovementioned methods are abbreviated as O-BEMD, S-BEMD, and T-BEMD, respectively. In this paper, the different BEMD algorithms are tested on several 256 256 standard images, including ‘‘LENA’’, ‘‘BARBARA’’, ‘‘RADISH’’, ‘‘PEPPERS’’, ‘‘TILE’’ and ‘‘FLOOR’’ (see the first column of Fig. 6, respectively), in which ‘‘LENA’’ and ‘‘BARBARA’’ are character images, ‘‘RADISH’’ and ‘‘PEPPERS’’ are plant images, whereas ‘‘TILE’’ and ‘‘FLOOR’’ are texture images. Additionally, we conduct qualitative and quantitative comparisons on those typical BEMD algorithms. 4.1. Boundary extension In this section, the experimental images are preprocessed by different boundary extension approaches according to the above-mentioned BEMD algorithms. More specifically, the O-BEMD does not take any measures to extend the boundary of the given image, and directly decomposes the image into several BIMFs. As mentioned before, the S-BEMD extends the given image symmetrically, whereas the T-BEMD involves a boundary extension technique inspired by nonparametric sampling based texture synthesis methods. Moreover, the proposed
D-BEMD implements Algorithm 1 of Section 3 for boundary extension. The processing order of each boundary is illustrated in Fig. 5. That is, both the S-BEMD and D-BEMD extend the boundary of each image in order of A-B-C-D (see Fig. 5(a)), whereas the T-BEMD synthesizes each pixel surrounding the original image from inside to outside (see Fig. 5(b)). As referred to [18], it is sufficient to extend 8 or 16 pixels (i.e. Q¼8 or 16) for most images. Particularly, we extend 16 pixels with T ¼ 15, e1 ¼ 0:1, e2 ¼ 0:15 for the six standard images based on different boundary extension methods utilized in S-BEMD, T-BEMD and D-BEMD. All the BEMD algorithms were tested in MATLAB R2010b, running on a Windows XP operating system with Intel Core2 CPU 2.93 GHz and 3.00 GB RAM. For the 256 256 images, the time complexity (i.e. the order of magnitude of the processing time) of boundary extension for different methods is displayed in Table 1. It is observed that S-BEMD is the fastest, D-BEMD is next, and T-BEMD is the slowest. This is due to the reason that the boundary processing algorithm in S-BEMD and DBEMD are patches based, whereas those in T-BEMD are pixels based, which are time-consuming. As depicted in Fig. 6, the first column shows the original images, which are the standard ‘‘LENA’’, ‘‘BARBARA’’, ‘‘RADISH’’, ‘‘PEPPERS’’, ‘‘TILE’’ and ‘‘FLOOR’’, respectively. The second to fourth columns illustrate decomposition results of S-BEMD, T-BEMD and D-BEMD, respectively. More precisely, 16 pixels closest to the border are symmetrically mapped to the external parts of the original image in SBEMD, which is a promising approach for symmetrical images (e.g. for the image ‘‘TILE’’). Nonetheless, the S-BEMD is not suitable for handling the asymmetrical case since many details of the symmetrical extension results (e.g. the mirror and hair in ‘‘LENA’’, the collar and arm in ‘‘BARBARA’’,
130
Z. He et al. / Signal Processing 93 (2013) 124–138
Begin Input IMoriginal(x,y) Boundary extension via DMGM(1,3), and obtain IMextended(x,y) Let r1=IMextended (x, y) ; i=1; j=0 j=j+1; hi(j-1)=ri
j=0;i=i+1 N
Identify maxima and minima of hi(j-1)
Y
Is ri+1 monotone?
Get envelopes emax , emin
Extract corresponding parts of BIMFextended
Compute mean Ave=(emax +emin )/2 ri+1=ri-ci
end hij=hi(j-1)-Ave
N
Y Is hij an BIMF?
Save the extended BIMFextended: ci=hij
Fig. 4. Flowchart of the proposed D-BEMD method. ‘‘A’’ and ‘‘B’’ demonstrate the major differences between D-BEMD and the traditional one.
D
IMoriginal C
IMoriginal B
A
Fig. 5. Boundary extension orders of different BEMD methods. (a) S-BEMD and D-BEMD: A-B-C-D. (b) T-BEMD: from inside to outside.
the underside leaves and fruit in ‘‘RADISH’’, the upper peppers and left peppers in ‘‘PEPPERS’’, as well as the texture in ‘‘FLOOR’’) cannot fully reflect the developing trend of the original images. Images in the third column are extended from nonparametric sampling based texture synthesis method. It is notable that this method does not always reveal excellent performance even though it verifies the trend of the original image to a certain extent.
For example, there are some unexpected complex textures on the left of Fig. 6(c), whose mirror is also extended unsatisfactorily. The complex details of Fig. 6(g) become smooth in the extended one, whereas the lower leaves blend into each other in Fig. 6(k). Moreover, Fig. 6(o) in the third column does not accurately represent the trend of the upper peppers. The detailed information is not described effectively in Fig. 6(s). Worse still, it appears confused
Z. He et al. / Signal Processing 93 (2013) 124–138
Table 1 Time complexity of boundary extension in various BEMD methods. Method
Time complexity (s)
S-BEMD
102
T-BEMDa
104
T-BEMDb
103
D-BEMD
101
a The original nonparametric sampling based method (referred to [19]). b The modified one (referred to [18]).
around Fig. 6(w). As illustrated in the fourth column of Fig. 6, the boundaries extended by Algorithm 1 significantly detect the developing trend of original images and yield considerable improvements compared with the former two methods.
131
Fig. 7(g). In fact, it turns out that the decomposition results of ‘‘LENA’’ from D-BEMD present better performance than other three methods. That is, the first BIMF (see Fig. 7(d)) has more detail information, and the second BIMF (see Fig. 7(h)) reflects the contour information more accurately. Similarly, Figs. Figs. 8–12 depict the results of different BEMD methods for ‘‘BARBARA’’, ‘‘RADISH’’, ‘‘PEPPERS’’, ‘‘TILE’’ and ‘‘FLOOR’’, respectively. According to the comparative study of the decomposition results, three major conclusions may be identified: (1) all the S-BEMD, T-BEMD and D-BEMD yield better results than that of the O-BEMD, which indicates that the boundary extension is necessary; (2) S-BEMD outperforms T-BEMD in case the image is symmetrical (e.g. ‘‘TILE’’), whereas the results gained from S-BEMD are similar to that of D-BEMD; (3) T-BEMD achieves better performance than S-BEMD in dealing with asymmetric image, whose results are less effective (i.e. ‘‘LENA’’, ‘‘BARBARA’’, ‘‘RADISH’’, ‘‘PEPPERS’’) or no less (i.e. ‘‘FLOOR’’) than that of D-BEMD. 4.3. Quantitative comparisons
4.2. Qualitative comparisons Having established the boundary extension of original images, the four different BEMD methods are applied to those modified images. As for earlier publications [4,12], two key issues among the traditional BEMD methods are stop criterion and interpolation strategy. In view of the stop criterion, we evaluate the standard deviation (SD) between two consecutive sifting results as
A further interesting issue addressed in this paper is to validate the proposed D-BEMD method quantitatively. 2 Given the original image IM original ðx,yÞ 2 l ðR2 Þ, ðx,y 2 þ N Þ with size P P, the decomposition result of a particular BEMD (such as the O-BEMD, S-BEMD, T-BEMD or D-BEMD) yields IM original ðx,yÞ ¼
I X
ci ðx,yÞ
ð15Þ
i¼1
SDij ¼
PX þ 2Q PX þ 2Q
2
9hði1Þj ðx,yÞhij ðx,yÞ9 2
x¼1 y¼1
hði1Þj ðx,yÞ
ð14Þ
where hij ðx,yÞ denotes the jth medium variable generated in finding the ith BIMF, whereas P þ 2Q gives the length (also the width) of the hij ðx,yÞ. On the other hand, followed by detecting the extrema via neighboring window methods (3 3 window is considered in this paper), Delaunay triangulation and cubic interpolation on triangles other than thin-plate splines are employed as interpolation scheme to decrease the computation time of BEMD. Consequently, the corresponding parts of the extended BIMFs together with the final residue are extracted as the eventual results of various methods. As displayed in Figs. 7–12, we extract the former two BIMFs and compare the decomposition results qualitatively. The first to fourth columns in Fig. 7 present the decomposition results of ‘‘LENA’’ obtained from O-BEMD, S-BEMD, T-BEMD and D-BEMD, respectively. It is observed from the first column that there exist serious boundary effects due to the absence of boundary effects mitigation measures in O-BEMD. More precisely, some vague pixels revealed on the left of Fig. 7(a), the second BIMF (see Fig. 7(e)) and residue (see Fig. 7(i)) represent obvious boundary effects. Indeed, it turns out that the other three BEMD methods establish improved results by virtue of various boundary effects alleviation strategies. Nevertheless, there remain ambiguous results in Fig. 7(f) obtained by S-BEMD. Moreover, undesired black blocks appear on the shoulder of the people in
where 1 rx, y rP, ci ðx,yÞ is the ith BIMF, I is supposed to be the total number of BIMFs, and the residue r I ðx,yÞ is taken as an additional element followed by [1]. Consequently, the square of the image IM original ðx,yÞ can be formed as !2 I X 2 ci ðx,yÞ IM original ðx,yÞ ¼ i¼1
¼
I X i¼1
c2i ðx,yÞ þ 2
I I X X
ci ðx,yÞcj ðx,yÞ
ð16Þ
i ¼ 1 j ¼ 1,jai
In case the decomposition result is ideal, the BIMFs will be orthogonal to each other and the cross terms from last part of Eq. (16) should be zero. Regarding these facts, we define the bi-dimensional index of orthogonality (BIO) as below: ! I I P X P X X X 9ci ðx,yÞcj ðx,yÞ9 BIO ¼ ð17Þ 2 i ¼ 1 j ¼ 1,jai x ¼ 1 y ¼ 1 IM original ðx,yÞ where 9 9 denotes the absolute, ci ðx,yÞcj ðx,yÞ gives the dot product between ci ðx,yÞ and cj ðx,yÞ, and IM2original ðx,yÞ presents the dot product of IM original ðx,yÞ and itself. In particular, the BIO reflects the overall orthogonality between BIMFs, which enables us to implement quantitative comparison of the above-mentioned BEMD methods. Moreover, it is worth noting that the nearer the BIO approximates to zero, the better the orthogonality of BIMFs. Table 2 demonstrates the BIO obtained from different BEMD methods. Since there are no boundary effects
132
Z. He et al. / Signal Processing 93 (2013) 124–138
Fig. 6. Boundary extension results of various BEMD methods for 256 256 standard images. The first column: original images, from top to bottom are ‘‘LENA’’, ‘‘BARBARA’’, ‘‘RADISH’’, ‘‘PEPPERS’’, ‘‘TILE’’ and ‘‘FLOOR’’. The second to fourth column: extended results of S-BEMD, T-BEMD and D-BEMD, respectively.
mitigation strategies in O-BEMD, evident boundary effects appear in the BIMFs and the border pixels tend to infinity, which leads to infinite BIO. That means, the BIO of
O-BEMD are invalid. Consequently, we compare the other three BEMD methods in aspect of the orthogonality. By analyzing Table 2, we have to remark that, the BIO
Z. He et al. / Signal Processing 93 (2013) 124–138
133
Fig. 7. Decomposition results of different BEMD methods for ‘‘LENA’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
Fig. 8. Decomposition results of different BEMD methods for ‘‘BARBARA’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
134
Z. He et al. / Signal Processing 93 (2013) 124–138
Fig. 9. Decomposition results of different BEMD methods for ‘‘RADISH’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
Fig. 10. Decomposition results of different BEMD methods for ‘‘PEPPERS’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
Z. He et al. / Signal Processing 93 (2013) 124–138
135
Fig. 11. Decomposition results of different BEMD methods for ‘‘TILE’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
Fig. 12. Decomposition results of different BEMD methods for ‘‘FLOOR’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
136
Z. He et al. / Signal Processing 93 (2013) 124–138
Appendix A. Proof of Theorem 1
Table 2 BIO of different BEMD methods. Standard images
LENA BARBARA RADISH PEPPERS TILE FLOOR
Different BEMD methods O-BEMD
S-EMD
T-BEMD
D-BEMD
1 1 1 1 1 1
0.0304 0.0116 0.0681 0.0361 0.0714 0.0576
0.0219 0.0159 0.0530 0.0227 0.0903 0.0532
0.0194 0.0110 0.0481 0.0211 0.0693 0.0455
In this appendix we indicate the proof of Theorem 1. To prove Theorem 1, we need first to discretize the PN1 ð1Þ ð1Þ C-GMC(1,N). Let cðtÞ ¼ dx1 ðtÞ=dt ¼ axð1Þ i ¼ 1 bi xi þ 1 1 ðtÞ þ ðtÞ þ bN , we derive Z k Z k ð1Þ dx1 ðtÞ cðtÞ dt ðA:1Þ dt ¼ dt k1 k1 which can be simplified to Z k ð1Þ ðkÞx ðk1Þ ¼ cðtÞ dt xð1Þ 1 1
ðA:2Þ
k1
obtained from T-BEMD is less than S-BEMD except the images are symmetric (e.g. ‘‘TILE’’) or complicated (e.g. ‘‘BARBARA’’), while D-BEMD deserves the smallest BIO in the three BEMD methods. As discussed above, the BIMFs of D-BEMD involves best orthogonality, while T-BEMD is next, and S-BEMD is worst. In a nutshell, the proposed D-BEMD achieves better performance than the other existing methods, including O-BEMD, S-BEMD, as well as T-BEMD.
Employing the trapezoid formula to Eq. (A.2) can be written explicitly as
5. Conclusions
xð1Þ 1 ðkÞ ¼
Aiming at mitigating boundary effects of BEMD, we have developed a novel boundary extension algorithm inspired by the discrete multivariate gray model. By theoretically proving the equivalence between the continuous and discrete multivariate gray model in Theorem 1, the D-GMC(1,3) is adopted to predict the inner pixels of each extended block in terms of existing border, whose relative data series and characteristic data series are obtained from the coordinates and pixels of the image, respectively. In comparison with existing techniques, including O-BEMD, S-BEMD, as well as T-BEMD, the contribution of our proposed D-BEMD can be outlined in two aspects. On one hand, the boundaries extended by Algorithm 1 accurately reflects the developing trend of original images and yields significant improvements compared with symmetrical extension or texture synthesis based techniques. On the other hand, as stated in Sections 4.2 and 4.3, D-BEMD establishes qualitative and quantitative improvements within reasonable time compared to other three BEMD methods mentioned above. Future research can be conducted, for instance, on the refinement of D-GMC(1,3) to improve the effectiveness of the D-BEMD, the establishment of the algorithmic evaluation system, the decomposition and de-noising of noise images, etc.
Acknowledgment The authors gratefully acknowledge the financial support of the National Basic Research Program of China (973 Program) (2012CB720000), the National Natural Science Foundation of China (Nos. 60975009 and 61171197), IBRSEM, HIT.KLOF.2010017 and Research Fund for the Doctoral Program of Higher Education of China (Nos. 20092302110037 and 20102302110033).
Rk
k1
cðtÞ dt,
1 x1ð1Þ ðkÞx1ð1Þ ðk1Þ ¼ ½cðkÞ þ cðk1Þ 2 " # 1 N 1 X 1 NX ð1Þ ð1Þ ð1Þ ¼ bi xið1Þ ðkÞ þ b x ðk1Þ þ 2b ax ðkÞax ðk1Þ N i iþ1 1 1 þ1 2 i¼1 i¼1
ðA:3Þ Consequently, we can deduce from Eq. (A.3) that N 1 X 10:5a ð1Þ bi x1 ðk1Þ þ xð1Þ i þ 1 ðk1Þ 1þ 0:5a 2ð1 þ 0:5aÞ i¼1
þ xð1Þ i þ 1 ðkÞÞþ
bN 1þ 0:5a
ðA:4Þ
It is interesting to see that Eq. (6) and Eq. (A.4) appear identical form, at this point, we assume the following expression: 8 10:5a > > < a1 ¼ 1 þ 0:5a ðA:5Þ bi > > : bi ¼ , i ¼ 1,2, . . . ,N 1þ 0:5a As previously mentioned, the predicted series of C-GMC(1,N) is evaluated as bN ð1eaðk1Þ Þ a Z k N 1 X þ eaðktÞ bi xð1Þ ðtÞ dt iþ1
ð1Þ aðk1Þ x^ 1 ðkÞ ¼ xð1Þ þ 1 ð1Þe
1
i¼1
aðk1Þ ¼ xð1Þ þ 1 ð1Þe
þ
N1 X i¼1
bi
Z
k 1
bN ð1eaðk1Þ Þ a
eaðktÞ xð1Þ i þ 1 ðtÞ dt,
k ¼ 2,3, . . . ,n þm ðA:6Þ
As referred to [21], the background value of the associated series xð1Þ ðtÞ, t 2 ½k1,k, k ¼ 2,3, . . . ,n þ m, iþ1 i ¼ 1,2, . . . ,N1 can be taken as 12 ðxð1Þ ðk1Þ þxð1Þ ðkÞÞ, iþ1 iþ1 According to Eq. (A.6), it turns out that ð1Þ aðk1Þ x^ 1 ðkÞ xð1Þ þ 1 ð1Þe
þ
N 1 X i¼1
bi
k2 X 1
2 j¼0
bN ð1eaðk1Þ Þ a
ðkjÞ þ xð1Þ ðkj1ÞÞ ðxð1Þ iþ1 iþ1
aðk1Þ ¼ xð1Þ þ 1 ð1Þe
Z
kj
eaðktÞ dt
kj1
N 1 k2 X X bN ð1ea Þeaj ð1eaðk1Þ Þ þ bi a 2a i¼1 j¼0
Z. He et al. / Signal Processing 93 (2013) 124–138
137
Table A1 The relationship between a and z. Values of a Values of z
0.1 0.000075
0.2 0.000533
0.3 0.001575
0.4 0.003200
0.5 0.005208
Values of a Values of z
0.6 0.007200
0.7 0.008575
0.8 0.008533
0.9 0.006075
1 0
ð1Þ ðxð1Þ i þ 1 ðkjÞ þxi þ 1 ðkj1ÞÞ,
k ¼ 2,3, . . . ,n þ m
ðA:7Þ
On the other hand, the predicted series of D-GMC(1,N) can be deduced from ð1Þ x^ 1 ðkÞ
¼ a
ð1Þ 1 x1 ðk1Þ þ
N 1 X
bi 2
i¼1
¼ a1 a1 xð1Þ 1 ðk2Þ þ
i¼1
þ
N1 X
bi
i¼1
2
bi 2
¼ a1 a1 ða
ð1Þ 1 x1 ðk3Þ þ
! ðk2Þ þxð1Þ ðk1ÞÞ þ bN ðxð1Þ iþ1 iþ1
N1 X i¼1
þ
N1 X
bi
i¼1
þ
N1 X
2
bi
i¼1
2
bi 2
ðk3Þ þxð1Þ ðk2ÞÞ þ bN Þ ðxð1Þ iþ1 iþ1
þ
þ
þ
ða1 Þj
ða1 Þj bN ,
bi 2
ðkjÞ þxð1Þ ðkj1ÞÞ ðxð1Þ iþ1 iþ1
k ¼ 2,3, . . . ,nþ m
ðA:8Þ
j¼0
The Maclaurin expansion of ea from Eq. (A.7) and a1 from Eq. (A.8), respectively, yields ea ¼ 1a þ
a2 a3 an þ þð1Þn þ oðan Þ 2! 3! n!
ðA:9Þ
10:5a 1þ 0:5a a a2 an ¼ 1a 1 þ 2 þ þ ð1Þn n þ oðan Þ 2 2 2
a1 ¼
¼ 1a þ
a2 a3 an þ 1 þ þ ð1Þn þ 1 n þ oðan þ 1 Þ 2 22 2
ðA:10Þ
With regard to the small value of a, we calculate the difference between a1 and ea approximately as
z ¼ ða1 Þea
a3 23
þ
a3 a4 a4 a3 þa4 ¼ þ 3! 23 4! 12
ðA:11Þ
from which we obtain the numerical relationship between a and z (see Table A1). As a consequence, there exists a1 ea , which leads to bi b 10:5a b b ¼ i ½1ða1 Þ i 1ea , ¼ i 1 bi ¼ 1 þ0:5a 1 þ 0:5a a a a i ¼ 1,2, . . . ,N
ðA:12Þ
Note that a1 e , bi ðbi =aÞð1e Þ, i ¼ 1,2, . . . ,N Pk2 aj and ð1ea Þ ¼ 1eaðk1Þ , Eq. (A.8) can be j¼0e a
ð1Þ ðxð1Þ i þ 1 ðkjÞ þ xi þ 1 ðkj1ÞÞ
k2 X bN
a
bi ð1Þ ð1ea Þeaj ðxð1Þ i þ 1 ðkjÞ þ xi þ 1 ðkj1ÞÞ 2a
eaj ð1ea Þ
aðk1Þ ¼ xð1Þ 1 ð1Þe
ðk1Þ þxð1Þ ðkÞÞ þ bN ðxð1Þ iþ1 iþ1
a
N 1 X k2 X
þ
bi ð1ea Þeaj ðxð1Þ ðkjÞ þ xð1Þ ðkj1ÞÞ iþ1 iþ1 2a i¼1j¼0
þ
bN ð1eaðk1Þ Þ a
¼ ða
k2 X
N 1 X k2 X
j¼0
N 1 X k2 X
2
ða1 Þj bN
i¼1j¼0
k1 ð1Þ x1 ð1Þ 1Þ
þ
k2 X
bi
j¼0
# ðk2Þ þxð1Þ ðk1ÞÞ þ bN ðxð1Þ iþ1 iþ1
i¼1j¼0
ða1 Þj
aðk1Þ xð1Þ 1 ð1Þe
¼
þ
N 1 X k2 X i¼1j¼0
ðk1Þ þxð1Þ ðkÞÞ þ bN ðxð1Þ iþ1 iþ1
"
ð1Þ x^ 1 ðkÞ ¼ ða1 Þk1 xð1Þ 1 ð1Þ
þ
ðk1Þ þxð1Þ ðkÞÞ þ bN ðxð1Þ iþ1 iþ1 N1 X
rewritten as
ðA:13Þ
Up till now, it is notable that Eq. (A.7) is the same as Eq. (A.13), which verifies the equivalence between C-GMC(1,N) and D-GMC(1,N). This completes the outline of proof of Theorem 1. References [1] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, E.H. Shih, Q. Zheng, C.C. Tung, H.H. Liu, The empirical mode decomposition method and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proceedings of the Royal Society of London, Series A (1998) 903–995. [2] Z.G. Xu, B.X. Huang, F. Zhang, Improvement of empirical mode decomposition under low sampling rate, Signal Processing 89 (11) (2009) 2296–2303. [3] N. Tsakalozos, K. Drakakis, S. Rickard, A formal study of the nonlinearity and consistency of the empirical mode decomposition, Signal Processing 92 (9) (2012) 1961–1969. [4] J.-C. Nunes, Y. Bouaoune, E. Delechelle, O. Niang, P. Bunel, Image analysis by bidimensional empirical mode decomposition, Image and Vision Computing 21 (12) (2003) 1019–1026. [5] S.M.A. Bhuiyan, N.O. Attoh-Okine, K.E. Barner, A.Y. Ayenu-Prah, R.R. Adhami, Bidimensional empirical mode decomposition using various interpolation techniques, Advances in Adaptive Data Analysis 1 (2) (2009) 309–338. [6] Z.H. Wu, N.E. Huang, X.Y. Chen, The multi-dimensional ensemble empirical mode decomposition method, Advances in Adaptive Data Analysis 1 (3) (2009) 339–372. [7] M.-H. Yeh, The complex bidimensional empirical mode decomposition, Signal Processing 92 (2) (2012) 523–541. [8] J. Fleureau, A. Kachenoura, L. Albera, J.-C. Nunes, L. Senhadji, Multivariate empirical mode decomposition and application to multichannel filtering, Signal Processing 91 (12) (2011) 2783–2792. [9] Y. Xu, B. Liu, J. Liu, S. Riemenschneider, Two-dimensional empirical mode decomposition by finite elements, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science (2006) 3081–3096.
138
Z. He et al. / Signal Processing 93 (2013) 124–138
[10] S.M.A. Bhuiyan, R.R. Adhami, J.F. Khan, Fast and adaptive bidimensional empirical mode decomposition using order-statistics filter based envelope estimation, EURASIP Journal on Advances in Signal Processing 2008 (21) (2008) 1–18. [11] R. Pabel, R. Koch, G. Jager, A. Kunoth, Fast empirical mode decompositions of multivariate data based on adaptive spline-wavelets and a generalization of the Hilbert–Huang-transform (HHT) to arbitrary space dimensions, Advances in Adaptive Data Analysis 2 (3) (2010) 337–358. [12] A. Linderhed, Adaptive image compression with wavelet packets and empirical mode decomposition, Department of Electrical Engi¨ neering, Linkoping University, 2004. [13] A. Linderhed, Image empirical mode decomposition: a new tool for image processing, Advances in Adaptive Data Analysis 1 (2) (2009) 265–294. [14] Y. Tian, K. Zhao, Y.P. Xu, F.Y. Peng, An image compression method based on the multi-resolution characteristics of BEMD, Computers & Mathematics with Applications 61 (8) (2011) 2142–2147. [15] P.-H. Chen, Y.-C. Yang, L.-M. Chang, Illumination adjustment for bridge coating images using BEMD-Morphology approach (BMA), Automation in Construction 19 (4) (2010) 475–484. [16] J.N. Huang, B.B. Zhao, Y.Q. Chen, P.D. Zhao, Bidimensional empirical mode decomposition (BEMD) for extraction of gravity anomalies associated with gold mineralization in the Tongshi gold field,
[17]
[18]
[19]
[20]
[21] [22] [23]
[24]
Western Shandong Uplifted Block, Eastern China, Computers & Geosciences 36 (7) (2010) 987–995. S.M.A. Bhuiyan, J.F. Khan, R.R. Adhami, A novel approach of edge detection via a fast and adaptive bidimensional empirical mode decomposition method, Advances in Adaptive Data Analysis 2 (2) (2010) 171–192. Z.X. Liu, S.L. Peng, Boundary processing of bidimensional EMD using texture synthesis, IEEE Signal Processing Letters 12 (1) (2005) 33–36. A.A. Efros, T.K. Leung, Texture synthesis by non-parametric sampling, in: Proceedings of the Seventh IEEE International Conference on Computer Vision, 1999, pp. 1033–1038. N.-M. Xie, S.-F. Liu, Discrete grey forecasting model and its optimization, Applied Mathematical Modelling 33 (2) (2009) 1173–1186. T.-L. Tien, A research on the grey prediction model GM(1,N), Applied Mathematics and Computation 218 (9) (2012) 4903–4916. T.-L. Tien, The indirect measurement of tensile strength by the new model FGMC (1,N), Measurement 44 (10) (2011) 1884–1897. Z. He, Y. Shen, Q. Wang, Boundary extension for Hilbert–Huang transform inspired by gray prediction model, Signal Processing 92 (3) (2012) 685–697. ¨ V. Kwatra, A. Schodl, I. Essa, G. Turk, A. Bobick, Graphcut textures: image and video synthesis using graph cuts, ACM Transactions on Graphics 22 (3) (2003) 277–286.
Contents lists available at SciVerse ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Discrete multivariate gray model based boundary extension for bi-dimensional empirical mode decomposition Zhi He, Qiang Wang, Yi Shen n, Yan Wang School of Astronautics, Harbin Institute of Technology (HIT), Harbin 150001, China
a r t i c l e in f o
abstract
Article history: Received 9 December 2011 Received in revised form 24 May 2012 Accepted 9 July 2012 Available online 20 July 2012
The bi-dimensional empirical mode decomposition (BEMD) has attracted extensive attention recently by virtue of its high performance in adaptive image processing. Unfortunately, this promising technique does not necessarily yield fruitful results due to the boundary effects. Motivated by the discrete multivariate gray model, we propose a boundary extension framework for mitigating the boundary effects of BEMD. In greater detail, followed by verifying the equivalence between the continuous and discrete multivariate gray model theoretically, a first-order three-variable discrete multivariate gray model D-GMC(1,3), which is derived from the continuous multivariate gray model with convolution integral C-GMC(1,N), is utilized to predict the middle pixels of each extended block in terms of existing border. Specifically, the coordinates and pixels of the image are respectively regarded as relative data series and characteristic data series of D-GMC(1,3). Experimental results on a set of widely used images indicate that the proposed approach can achieve qualitative and quantitative improvements within appropriate processing time by comparing with other three generally acknowledged methods, i.e. the original BEMD, symmetrical extension as well as texture synthesis based BEMD. & 2012 Elsevier B.V. All rights reserved.
Keywords: Bi-dimensional empirical mode decomposition (BEMD) Boundary effects Discrete multivariate gray model
1. Introduction Over the last decades, the bi-dimensional empirical mode decomposition (BEMD) has become one of the most significant image processing techniques owing to its adaptive nature and versatile capability. In analogy to the one-dimensional (1-D) empirical mode decomposition (EMD) [1–3], the BEMD is also adaptive, empirical and intuitive due to its fully data-driven property. To date, the development of BEMD is twofold: theory and application. More precisely, Nunes et al. developed the EMD into two-dimensional (2-D) version, namely BEMD [4], in which the morphological reconstruction and radial basis
n Corresponding author. Tel.: þ 86 451 86413411; fax: þ86 451 86418378. E-mail address: [email protected] (Y. Shen).
0165-1684/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sigpro.2012.07.009
function (RBF) are implemented to detect the image extrema and calculate the surface interpolation, respectively. Bhuiyan et al. evaluated the existing surface interpolation techniques and concluded that the RBF based interpolators show promising performance [5]. To alleviate the scale mixing, a noise-assisted method termed as multi-dimensional ensemble EMD (MEEMD) was introduced by Wu et al. [6]. Recently, Yeh established a complex extension of BEMD, which took advantage of the perfect reconstruction of the four quadrant spectra and suitable for 2-D complex signal processing [7]. Fleureau et al. studied the multivariate EMD in terms of redefining the local mean operator [8]. Moreover, a number of fast algorithms based on finite elements [9], order-statistics filter [10] or adaptive spline-wavelets [11] have also been exploited. On the other hand, different kinds of applications are booming up nowadays. For instance, Linderhed [12,13] and Tian et al. [14] developed
Z. He et al. / Signal Processing 93 (2013) 124–138
a totally different image compression approach inspired by the multi-resolution properties of bi-dimensional intrinsic mode function (BIMF). A BEMD-morphology approach (BMA), which was proposed by Chen et al., is effective in handling non-uniform illumination in bridge coating assessment [15]. Huang et al. applied the BEMD to decompose the gravity data for the Tongshi gold field and obtained a clearer spatial distribution than that of Fourier transform [16]. As stated in [17], the desired edges of image were detected by a fast and adaptive BEMD algorithm abbreviated as FABEMD. Though BEMD has recently developed a valuable tool for image processing, how to mitigate boundary effects (see Fig. 1) remains an open problem, generated by the lack of extrema on the boundary or beyond the original range of the image. Indeed, the extrema may be sparse thus the envelopes detected by interpolating local extrema cannot cover the existing region of the image. Worse still, those drawbacks may propagate into the interior part of the image and ‘‘pollute’’ the full decomposition results. Much work has been carried out in the literature to handle this issue, whose merit is to extend data explicitly or implicitly out of the original region as suggested by Huang et al. in 1-D case [1]. Typically, the
125
symmetrical or periodic extension based BEMD extends the given image symmetrically or periodically. Nevertheless, the performance of this method is far from satisfactory in case the image is asymmetrical or nonperiodic. Moreover, Liu et al. introduced an extension technique inspired by nonparametric sampling based texture synthesis method [18,19]. Unfortunately, two main problems related to this approach: (1) it requires that the given image is extended accurately enough before decomposition; (2) it desires quantitative judgment of various forms of boundary processing. Other than the above-mentioned methods, we propose an alternative boundary extension method motivated by the discrete multivariate gray model [20–22]. Indeed, as a valuable forecasting technique, the gray model has successfully alleviated the end effects of Hilbert–Huang transform in our previous work [23]. More critically, in this paper, followed by evaluating the equivalence between the continuous and discrete multivariate gray model theoretically, a first-order three-variable discrete multivariate gray model termed as D-GMC(1,3), which is derived from the continuous multivariate gray model with convolution integral C-GMC(1,N), is associated with the extension of the given image. In general, we regard the coordinates and pixels of the image as relative data series and characteristic data series of D-GMC(1,3), respectively. Once the boundary of the given image is extended by D-GMC(1,3), BEMD is devoted to extracting the essential components, namely BIMFs. It is noteworthy that the bi-dimensional index of orthogonality (BIO) is introduced to evaluate various boundary processing methods quantitatively. The outline of the paper is organized as follows: Section 2 introduces the principle of the multivariate gray model. This is followed by a detailed description of our proposed boundary extension scheme in Section 3. In Section 4, comparative analysis of the widespread methods is carried out in terms of extensive experiments on standard images. The conclusion is drawn in Section 5. 2. The multivariate gray model 2.1. The continuous multivariate gray model
Fig. 1. An example of boundary effects during extracting the first BIMF. (a) The 3-D graphic of ‘‘LENA’’ depicted in Fig. 6(a). (b) The mean envelope in the process of obtaining the first BIMF.
As mentioned earlier, the boundary effects of BEMD are mainly due to lacking of extrema on the boundary or beyond the original region of the image. To alleviate this annoying deficiency, a boundary extension algorithm inspired by the discrete multivariate gray mode is proposed in this paper. Alternatively, the keynote of this section focuses on the overview of continuous multivariate gray model with convolution integral, termed as C-GMC(1,N) [21], which is recently proposed by Tien in view of the incorrect solution of the original gray model GM(1,N). Indeed, it is the foundation of discrete multivariate gray model as will be described in the following section. Assume the original characteristic data series and ð0Þ ð0Þ ð0Þ relative data series are Xð0Þ 1 ¼ fx1 ð1Þ,x1 ð2Þ, . . . ,x1 ðnÞg ð0Þ ð0Þ ð0Þ ð0Þ and Xi ¼ fxi ð1Þ,xi ð2Þ, . . . ,xi ðn þmÞg, i ¼ 2,3, . . . ,N, respectively, where m is the number of entries to
126
Z. He et al. / Signal Processing 93 (2013) 124–138
be predicted, xð0Þ and xð0Þ ðkÞ Z0, 1 ðkÞ Z 0,k ¼ 1,2, . . . ,n i ð0Þ k ¼ 1,2, . . . ,n þ m, i ¼ 2,3, . . . ,N. xi ðkÞ can be considered as integer values sampled from a continuous function. Pk ð0Þ ð1Þ Given xð1Þ l ¼ 1 x1 ðlÞ, k ¼ 1,2, . . . ,n and xi ðkÞ ¼ 1 ðkÞ ¼ Pk ð0Þ k ¼ 1,2, . . . ,n þ m, i ¼ 2,3, . . . ,N as the l ¼ 1 xi ðlÞ, first-order accumulated generating operation (1-AGO) of xið0Þ ðkÞ, the whitenization form of C-GMC(1,N) can then be defined as ð1Þ N 1 X dx1 ðtÞ bi xð1Þ ðtÞ þbN , þ axð1Þ 1 ðtÞ ¼ iþ1 dt i¼1
2.2. The discrete multivariate gray model Versatile as C-GMC(1,N), the improvement of prediction accuracy remains a challenging problem. So far, much work has been carried out in the literature to handle this issue [21,22]. Despite the partial improvements, all these methods do not solve the problem intrinsically. As stated in [20], the coefficients are evaluated by discrete equation, whereas the modeling and predicting are given by continuous equations. Consequently, the paradox between discrete estimation and continuous prediction is just the cause of forecast error in C-GMC(1,N). To counter this effect, we introduce the discrete multivariate gray model, in which the parameter estimation, modeling construction together with data predication are achieved by discrete equations and there is no contradiction between discrete estimation and continuous prediction, and therefore, has high prediction accuracy.
t ¼ 1,2, . . . ,n þ m ð1Þ
where a is gray developmental coefficient, bi ,i ¼ 1,2, . . . ,N are associated coefficients and rp ¼0 in Eq. (33) of [21] taking into account the actual situation of this paper. Letting k1 ¼ ½a b1 b2 . . . bN T , we can determine k1 as follows
k1 ¼ ½a b1 b2 . . . bN T ¼ ðBT1 B1 Þ1 BT1 Y 1
Definition 1. Assume the original characteristic data ð0Þ series Xð0Þ 1 and relative data series Xi , i ¼ 2,3, . . . ,N with n entries are the same as those of Section 2.1. The
ð2Þ
ð0Þ T where Y 1 ¼ ½x1ð0Þ ð2Þ xð0Þ 1 ð3Þ . . . x1 ðnÞ and
2
ð1Þ 12 ðxð1Þ 1 ð1Þ þ x1 ð2ÞÞ 6 6 1 ðxð1Þ ð2Þ þ xð1Þ ð3ÞÞ 6 1 2 1 B1 ¼ 6 6 ^ 4 ð1Þ 12 ðxð1Þ 1 ðn1Þ þx1 ðnÞÞ
ð1Þ 1 ð1Þ 2 ðx2 ð1Þ þ x2 ð2ÞÞ ð1Þ 1 ð1Þ 2 ðx2 ð2Þ þ x2 ð3ÞÞ
^
&
^
ð1Þ 1 ð1Þ 2 ðx2 ðn1Þ þx2 ðnÞÞ
ð1Þ 1 ð1Þ 2 ðxN ðn1Þ þxN ðnÞÞ
aðt1Þ ¼ xð1Þ þ 1 ð1Þe
aðt1Þ ¼ xð1Þ þ 1 ð1Þe
Z
Z
t
eaðttÞ
1 t
N 1 X
! bi xð1Þ ðtÞ þ bN iþ1
ð1Þ xð1Þ 1 ðkÞ ¼ a1 x1 ðk1Þ þ
dt
Z
1
t
eaðttÞ
1
bN ð1eaðt1Þ Þþ a
Z
t
eaðttÞ
1
N 1 X i¼1
N 1 X
bi xð1Þ ðtÞ dt iþ1
i¼1
ð3Þ
þ 1
aðktÞ
e
bN ð1eaðk1Þ Þ a
N 1 X
bi xð1Þ i þ 1ð
tÞ dt
ð4Þ
i¼1
ð1Þ
where x^ 1 ðkÞ denotes the predicted value of xð1Þ 1 ðkÞ. Therefore, the predicted value of x1ð0Þ ðkÞ, k ¼ 2,3, . . . ,n þ m yields ð0Þ ð1Þ ð1Þ x^ 1 ðkÞ ¼ x^ 1 ðkÞx^ 1 ðk1Þ
bi 2
ðxð1Þ ðk1Þ þ xð1Þ ðkÞÞ þ bN iþ1 iþ1 ð6Þ
bi xð1Þ ðtÞ dt iþ1
P ð1Þ where 1 eaðttÞ ð N1 i ¼ 1 bi xi þ 1 ðtÞ þ bN Þ dt is just the convolution integral that the C-GMC(1,N) differs from the original GM(1,N). Consequently, the predicted series of xð1Þ 1 ðkÞ,k ¼ 2,3, . . . ,n þ m are deduced from
k
N 1 X i¼1
Rt
Z
3
7 17 7 7 ^7 5 1
i¼1
eaðttÞ bN dt þ
ð1Þ aðk1Þ þ x^ 1 ðkÞ ¼ xð1Þ 1 ð1Þe
1
first-order N-variable discrete multivariate gray model, namely D-GMC(1,N), can be defined as
Solving Eq. (1) leads to ð1Þ aðt1Þ xð1Þ þ 1 ðtÞ ¼ x1 ð1Þe
ð1Þ 1 ð1Þ 2 ðxN ð1Þ þ xN ð2ÞÞ ð1Þ 1 ð1Þ 2 ðxN ð2Þ þ xN ð3ÞÞ
ð5Þ
Pk
ð0Þ where xð1Þ l ¼ 1 xi ðlÞ, k ¼ 1,2, . . . ,n þ m, i ¼ 1,2, . . . ,N i ðkÞ ¼ and k2 ¼ ½a1 b1 b2 . . . bN T is the coefficient.
In analogy to C-GMC(1,N), the coefficient k2 is evaluated as
k2 ¼ ½a1 b1 b2 . . . bN T ¼ ðBT2 B2 Þ1 BT2 Y 2 where
Y 2 ¼ ½xð1Þ 1 ð2Þ
2
xð1Þ 1 ð1Þ 6 6 xð1Þ ð2Þ 6 1 B2 ¼ 6 6 ^ 4 ð1Þ x1 ðn1Þ
xð1Þ 1 ð3Þ
...
T xð1Þ 1 ðnÞ
ð1Þ 1 ð1Þ 2 ðx2 ð1Þ þ x2 ð2ÞÞ ð1Þ 1 ð1Þ ðx ð2Þ þ x 2 ð3ÞÞ 2 2
^ ð1Þ 1 ð1Þ 2 ðx2 ðn1Þþ x2 ðnÞÞ
and ð1Þ 1 ð1Þ 2 ðxN ð1Þ þ xN ð2ÞÞ ð1Þ 1 ð1Þ ðx ð2Þ þ x N ð3ÞÞ 2 N
&
ð7Þ
^ ð1Þ 1 ð1Þ 2 ðxN ðn1Þ þ xN ðnÞÞ
1
3
7 17 7 7 ^7 5
1
Given the numerous coefficients of D-GMC(1,N), the following intermediate matrix U is introduced to avoid excessive matrix computations: 2
u11 6 u 6 21 U¼6 6 ^ 4 uN þ 1,1
u12
u22
^
&
uN þ 1,2
u1,N þ 1
3
7 u2,N þ 1 7 7 ¼ ðBT B2 Þ1 2 7 ^ 5
uN þ 1,N þ 1
ð8Þ
Z. He et al. / Signal Processing 93 (2013) 124–138
Substituting Eq. (8) into Eq. (7) leads to 2 3 2 3 2 a1 xð1Þ u12 u1,N þ 1 u11 1 ð1Þ 6b 7 6 17 6 7 6 6 u22 u2,N þ 1 7 6 xð1Þ 6 7 6 u21 1 ð2Þ 76 b 7 6 k2 ¼ 6 6 2 7¼6 ^ 7 6 ^ & ^ ^ 6 ^ 7 4 5 4 4 5 uN þ 1,1 uN þ 1,2 uN þ 1,N þ 1 xð1Þ ðn1Þ
bN
1
2
127
ð1Þ 1 ð1Þ 2 ðx2 ð1Þ þ x2 ð2ÞÞ ð1Þ 1 ð1Þ 2 ðx2 ð2Þ þ x2 ð3ÞÞ
^
&
^
ð1Þ 1 ð1Þ 2 ðx2 ðn1Þ þ x2 ðnÞÞ
ð1Þ 1 ð1Þ 2 ðxN ðn1Þ þxN ðnÞÞ
ð1Þ 1 ð1Þ 2 ðxN ð1Þ þ xN ð2ÞÞ ð1Þ 1 ð1Þ 2 ðxN ð2Þ þ xN ð3ÞÞ
3 n X ðxð1Þ ðk1Þ þxð1Þ ðkÞÞxð1Þ ð1Þ 1 ðkÞ i i þu1,N þ 1 u11 xð1Þ u1i xð1Þ 6 7 1 ðk1Þx1 ðkÞ þ 1 ðkÞ 2 7 2 ð1Þ 3 6 i¼2 k¼2 k¼2 k¼2 6 7 x1 ð2Þ 6 7 ð1Þ ð1Þ ð1Þ n N n n X X ð1Þ X X 6 ð1Þ 7 6 7 ðx ðk1Þ þx ðkÞÞx ðkÞ 1 i i 6 x ð3Þ 7 6 7 þu2,N þ 1 x1 ðk1Þxð1Þ u2i xð1Þ u21 6 1 7 6 7 1 ðkÞ þ 1 ðkÞ 2 6 7¼6 7 i¼2 k¼2 k¼2 k¼2 6 ^ 7 6 7 4 5 6 7 ^ 6 7 ð1Þ x1 ðnÞ 6 7 ð1Þ ð1Þ ð1Þ n N n n X X X X 6 7 ðxi ðk1Þ þxi ðkÞÞx1 ðkÞ ð1Þ ð1Þ ð1Þ 4 uN þ 1,1 þuN þ 1,N þ 1 x1 ðk1Þx1 ðkÞ þ uN þ 1,i x1 ðkÞ 5 2 i¼2 k¼2 k¼2 k¼2 n X
N X
Setting 8 n X > ð1Þ > > xð1Þ > CðjÞ ¼ uj1 1 ðk1Þx1 ðkÞ > > > k¼2 > > > > N n < X X ðxð1Þ ðk1Þ þxð1Þ ðkÞÞxð1Þ 1 ðkÞ i i uji DðjÞ ¼ 2 > > i ¼ 2 k ¼ 2 > > > n > X > > > EðjÞ ¼ u xð1Þ > j,N þ 1 1 ðkÞ > : k¼2 the coefficient k2 in Eq. (9) can be rewritten as 2 3 2 3 a1 Cð1Þ þ Dð1Þ þEð1Þ 6b 7 6 17 6 7 Cð2Þ þ Dð2Þ þEð2Þ 6 7 6 7 7¼6 7 b k2 ¼ 6 2 6 7 6 7 ^ 6 ^ 7 4 5 4 5 CðN þ 1Þ þ DðN þ 1Þ þEðN þ 1Þ
n X
ð9Þ
ð10Þ Algorithm 1. Boundary extension algorithm via D-GMC(1,3). 2
ð1Þ x^ ðkÞ ¼ a1 xð1Þ 1 ðk1Þ þ
N 1 X i¼1
bi 2
ðxð1Þ ðk1Þ þ xð1Þ ðkÞÞ þ bN iþ1 iþ1 ð11Þ
and ð0Þ ð1Þ ð1Þ x^ 1 ðkÞ ¼ x^ 1 ðkÞx^ 1 ðk1Þ
ð12Þ
ð0Þ and x^ 1 ðkÞ, k ¼ 2,3, . . . ,n þ m are ð1Þ of x1 ðkÞ and x1ð0Þ ðkÞ, respectively.
where the predicted values We can assess the forecast performance of C-GMC(1,N) or D-GMC(1,N) via the following average relative error: n 1X L nk¼1 k
7 17 7 7 ^7 5 1
Based on the study of the discrete multivariate gray model, we can summarize the proposed boundary extension algorithm in Algorithm 1. In greater detail, 2 the original image IM original ðx,yÞ 2 l ðR2 Þ, ðx,y 2 N þ Þ, whose size is P P, can be predicted by the D-GMC (1,3), in which the coordinates and pixels of the image are respectively regarded as relative data series and characteristic data series. As such, the size of extended image IM extended ðx,yÞ becomes ðP þ 2Q Þ ðP þ2Q Þ with the central P P block being the original image IM original ðx,yÞ.
Having decided the coefficient k2 , the predicted values ð0Þ of xð1Þ 1 ðkÞ and x1 ðkÞ,k ¼ 2,3, . . . ,n þ m can be determined by
L¼
3T
3. The proposed boundary extension algorithm
bN
ð1Þ x^ 1 ðkÞ
1
ð13Þ
where Lk ¼ 9eðkÞ=x1ð0Þ ðkÞ9, and eðkÞ ¼ x1ð0Þ ðkÞx^ ð0Þ 1 ðkÞ, k ¼ 1,2, . . . ,n. Finally, Theorem 1 verifies the the equivalence between C-GMC(1,N) and D-GMC(1,N). Theorem 1. Assume the gray developmental coefficient a of Eq. (1) is small, then C-GMC(1,N) is equivalent to D-GMC(1,N). Proof. The proof of Theorem 1 is given in Appendix A.
Require: The original image IM original ðx,yÞ 2 l ðR2 Þ,ðx,y 2 N þ Þ with size P P; error e1 and e2 ; size of the extended image ðP þ 2Q Þ ðP þ 2Q Þ; size of searching block T T Ensure The extended image IM extended ðx,yÞ 1: %initialization 2: IMprocess ðx,yÞ ¼ IMoriginal ðx,yÞ 3: %extension procedure 4: while IMprocess ðx,yÞ is not fully filled into ðP þ 2Q Þ ðP þ 2Q Þ do 5: Browse through the original image IMoriginal ðx,yÞ in raster scan order 6: Search IM original ðx,yÞ for a set of blocks (sizes are T T) B ¼ fB1 ,B2 , . . . ,Bl g that fulfill the overlap constraints from above and left (only the about overlap is considered in case filling the first column of IM process ðx,yÞ) within the error tolerance e1 7: Pick one block Brandom randomly from B ¼ fB1 ,B2 , . . . ,Bl g, B0 ¼ Brandom 8: Predict the inner pixel values of B0 via D-GMC(1,3) in terms of the surrounding pixel values of Brandom 9: Get the prediction error L of D-GMC(1,3) 10: 11: 12: 13: 14: 15: 16: 17:
if L 4 e2 then B0 ¼ Brandom %in case the prediction accuracy of D-GMC(1,3) is unsatisfactory end if Calculate the error surface S between the newly generated block B0 and the old blocks on the overlap region Identify the minimum cost path along S and set it as the boundary of B0 [24] Tile the modified block B0 onto the corresponding region of IMprocess ðx,yÞ end while IMextended ðx,yÞ ¼ IM process ðx,yÞ
128
Z. He et al. / Signal Processing 93 (2013) 124–138
characteristic data series X ð0Þ 1 and relative data series ð0Þ X ð0Þ and X of D-GMC(1,3) are given by pixel values, 2 3 horizontal as well as vertical coordinates, respectively, of the chosen block B0 . As long as the prediction error is smaller than e2 , these predicted inner pixels will be considered as a partial extended boundary of IM extended ðx,yÞ. 4. It is highlighted from Fig. 3 that the boundary of B0 is evaluated by the minimum cost path [24] along the error surface S. In this way, we can improve the smoothness and cohesive property of overlapping between two small blocks. In greater detail, Fig. 3 outlined all the three potential overlapping scenarios, horizontal, vertical and corner overlaps, respectively.
The following four remarks further clarify Algorithm 1: 1. As depicted in Fig. 2(a), we divide the boundary to be extended into eight small blocks. It is worth underlining that different combinations of those blocks may yield different extension orders, which are typically shown in Fig. 2(b)–(g). Although various combinations of the blocks may influence the boundary extension results, these influences are quite small in general. From this viewpoint, we choose the extension order according to Fig. 2(b) without loss of generality. 2. It is illustrated in Fig. 3 that the extension order of original image is A-B-C-D, which is combined from the eight small blocks of Fig. 2(a). As for ‘‘A’’, it is sequently filled by a few small blocks with sizes T T (see Fig. 3), each of them are randomly picked from a set of candidate blocks. Particularly, the other boundaries, i.e., ‘‘B’’, ‘‘C’’ and ‘‘D’’, can be converted into the same situation as ‘‘A’’ by virtue of folding, rotation as well as transposition. 3. With regard to each small block B0 in Fig. 3, its inner pixels can be predicted by D-GMC(1,3) making full use of all their surrounding pixels. More specifically, the
1
4
6
Having finished illustrating the proposed boundary extension approach, the modified BEMD algorithm inspired by D-GMC(1,3), namely D-BEMD, can be drawn in Fig. 4. It is noteworthy that the main contribution of D-BEMD is twofold:
The original image is extended by D-GMC(1,3) before applying the traditional BEMD.
2
IMoriginal
7
3
5
8
IMoriginal
IMoriginal
IMoriginal
IMoriginal
IMoriginal
IMoriginal
Fig. 2. Schematic division of the extended blocks in Algorithm 1. (a) Eight small blocks to be extended. (b)–(g) Six different combinations of the small blocks in (a).
Z. He et al. / Signal Processing 93 (2013) 124–138
129
Pick a block randomly
B1 B2 Bl
B0
IMoriginal
Enlarged
DMGM (1,3)
D Minimum cost path
A IMoriginal C
B Enlarged
Potential overlapping scenarios
A Horizontal Vertical
Corner
Fig. 3. Schematic presentation of Algorithm 1.
The corresponding parts of the extended BIMFs and final residue obtained from the extended image are extracted as decomposition results of the original image.
4. Experimental results The keynote of this section focuses on assessing the performance of the proposed D-BEMD algorithm by comparing with other three generally acknowledged methods, i.e. the original BEMD, symmetrical extension as well as texture synthesis based BEMD. For the sake of simplicity, the abovementioned methods are abbreviated as O-BEMD, S-BEMD, and T-BEMD, respectively. In this paper, the different BEMD algorithms are tested on several 256 256 standard images, including ‘‘LENA’’, ‘‘BARBARA’’, ‘‘RADISH’’, ‘‘PEPPERS’’, ‘‘TILE’’ and ‘‘FLOOR’’ (see the first column of Fig. 6, respectively), in which ‘‘LENA’’ and ‘‘BARBARA’’ are character images, ‘‘RADISH’’ and ‘‘PEPPERS’’ are plant images, whereas ‘‘TILE’’ and ‘‘FLOOR’’ are texture images. Additionally, we conduct qualitative and quantitative comparisons on those typical BEMD algorithms. 4.1. Boundary extension In this section, the experimental images are preprocessed by different boundary extension approaches according to the above-mentioned BEMD algorithms. More specifically, the O-BEMD does not take any measures to extend the boundary of the given image, and directly decomposes the image into several BIMFs. As mentioned before, the S-BEMD extends the given image symmetrically, whereas the T-BEMD involves a boundary extension technique inspired by nonparametric sampling based texture synthesis methods. Moreover, the proposed
D-BEMD implements Algorithm 1 of Section 3 for boundary extension. The processing order of each boundary is illustrated in Fig. 5. That is, both the S-BEMD and D-BEMD extend the boundary of each image in order of A-B-C-D (see Fig. 5(a)), whereas the T-BEMD synthesizes each pixel surrounding the original image from inside to outside (see Fig. 5(b)). As referred to [18], it is sufficient to extend 8 or 16 pixels (i.e. Q¼8 or 16) for most images. Particularly, we extend 16 pixels with T ¼ 15, e1 ¼ 0:1, e2 ¼ 0:15 for the six standard images based on different boundary extension methods utilized in S-BEMD, T-BEMD and D-BEMD. All the BEMD algorithms were tested in MATLAB R2010b, running on a Windows XP operating system with Intel Core2 CPU 2.93 GHz and 3.00 GB RAM. For the 256 256 images, the time complexity (i.e. the order of magnitude of the processing time) of boundary extension for different methods is displayed in Table 1. It is observed that S-BEMD is the fastest, D-BEMD is next, and T-BEMD is the slowest. This is due to the reason that the boundary processing algorithm in S-BEMD and DBEMD are patches based, whereas those in T-BEMD are pixels based, which are time-consuming. As depicted in Fig. 6, the first column shows the original images, which are the standard ‘‘LENA’’, ‘‘BARBARA’’, ‘‘RADISH’’, ‘‘PEPPERS’’, ‘‘TILE’’ and ‘‘FLOOR’’, respectively. The second to fourth columns illustrate decomposition results of S-BEMD, T-BEMD and D-BEMD, respectively. More precisely, 16 pixels closest to the border are symmetrically mapped to the external parts of the original image in SBEMD, which is a promising approach for symmetrical images (e.g. for the image ‘‘TILE’’). Nonetheless, the S-BEMD is not suitable for handling the asymmetrical case since many details of the symmetrical extension results (e.g. the mirror and hair in ‘‘LENA’’, the collar and arm in ‘‘BARBARA’’,
130
Z. He et al. / Signal Processing 93 (2013) 124–138
Begin Input IMoriginal(x,y) Boundary extension via DMGM(1,3), and obtain IMextended(x,y) Let r1=IMextended (x, y) ; i=1; j=0 j=j+1; hi(j-1)=ri
j=0;i=i+1 N
Identify maxima and minima of hi(j-1)
Y
Is ri+1 monotone?
Get envelopes emax , emin
Extract corresponding parts of BIMFextended
Compute mean Ave=(emax +emin )/2 ri+1=ri-ci
end hij=hi(j-1)-Ave
N
Y Is hij an BIMF?
Save the extended BIMFextended: ci=hij
Fig. 4. Flowchart of the proposed D-BEMD method. ‘‘A’’ and ‘‘B’’ demonstrate the major differences between D-BEMD and the traditional one.
D
IMoriginal C
IMoriginal B
A
Fig. 5. Boundary extension orders of different BEMD methods. (a) S-BEMD and D-BEMD: A-B-C-D. (b) T-BEMD: from inside to outside.
the underside leaves and fruit in ‘‘RADISH’’, the upper peppers and left peppers in ‘‘PEPPERS’’, as well as the texture in ‘‘FLOOR’’) cannot fully reflect the developing trend of the original images. Images in the third column are extended from nonparametric sampling based texture synthesis method. It is notable that this method does not always reveal excellent performance even though it verifies the trend of the original image to a certain extent.
For example, there are some unexpected complex textures on the left of Fig. 6(c), whose mirror is also extended unsatisfactorily. The complex details of Fig. 6(g) become smooth in the extended one, whereas the lower leaves blend into each other in Fig. 6(k). Moreover, Fig. 6(o) in the third column does not accurately represent the trend of the upper peppers. The detailed information is not described effectively in Fig. 6(s). Worse still, it appears confused
Z. He et al. / Signal Processing 93 (2013) 124–138
Table 1 Time complexity of boundary extension in various BEMD methods. Method
Time complexity (s)
S-BEMD
102
T-BEMDa
104
T-BEMDb
103
D-BEMD
101
a The original nonparametric sampling based method (referred to [19]). b The modified one (referred to [18]).
around Fig. 6(w). As illustrated in the fourth column of Fig. 6, the boundaries extended by Algorithm 1 significantly detect the developing trend of original images and yield considerable improvements compared with the former two methods.
131
Fig. 7(g). In fact, it turns out that the decomposition results of ‘‘LENA’’ from D-BEMD present better performance than other three methods. That is, the first BIMF (see Fig. 7(d)) has more detail information, and the second BIMF (see Fig. 7(h)) reflects the contour information more accurately. Similarly, Figs. Figs. 8–12 depict the results of different BEMD methods for ‘‘BARBARA’’, ‘‘RADISH’’, ‘‘PEPPERS’’, ‘‘TILE’’ and ‘‘FLOOR’’, respectively. According to the comparative study of the decomposition results, three major conclusions may be identified: (1) all the S-BEMD, T-BEMD and D-BEMD yield better results than that of the O-BEMD, which indicates that the boundary extension is necessary; (2) S-BEMD outperforms T-BEMD in case the image is symmetrical (e.g. ‘‘TILE’’), whereas the results gained from S-BEMD are similar to that of D-BEMD; (3) T-BEMD achieves better performance than S-BEMD in dealing with asymmetric image, whose results are less effective (i.e. ‘‘LENA’’, ‘‘BARBARA’’, ‘‘RADISH’’, ‘‘PEPPERS’’) or no less (i.e. ‘‘FLOOR’’) than that of D-BEMD. 4.3. Quantitative comparisons
4.2. Qualitative comparisons Having established the boundary extension of original images, the four different BEMD methods are applied to those modified images. As for earlier publications [4,12], two key issues among the traditional BEMD methods are stop criterion and interpolation strategy. In view of the stop criterion, we evaluate the standard deviation (SD) between two consecutive sifting results as
A further interesting issue addressed in this paper is to validate the proposed D-BEMD method quantitatively. 2 Given the original image IM original ðx,yÞ 2 l ðR2 Þ, ðx,y 2 þ N Þ with size P P, the decomposition result of a particular BEMD (such as the O-BEMD, S-BEMD, T-BEMD or D-BEMD) yields IM original ðx,yÞ ¼
I X
ci ðx,yÞ
ð15Þ
i¼1
SDij ¼
PX þ 2Q PX þ 2Q
2
9hði1Þj ðx,yÞhij ðx,yÞ9 2
x¼1 y¼1
hði1Þj ðx,yÞ
ð14Þ
where hij ðx,yÞ denotes the jth medium variable generated in finding the ith BIMF, whereas P þ 2Q gives the length (also the width) of the hij ðx,yÞ. On the other hand, followed by detecting the extrema via neighboring window methods (3 3 window is considered in this paper), Delaunay triangulation and cubic interpolation on triangles other than thin-plate splines are employed as interpolation scheme to decrease the computation time of BEMD. Consequently, the corresponding parts of the extended BIMFs together with the final residue are extracted as the eventual results of various methods. As displayed in Figs. 7–12, we extract the former two BIMFs and compare the decomposition results qualitatively. The first to fourth columns in Fig. 7 present the decomposition results of ‘‘LENA’’ obtained from O-BEMD, S-BEMD, T-BEMD and D-BEMD, respectively. It is observed from the first column that there exist serious boundary effects due to the absence of boundary effects mitigation measures in O-BEMD. More precisely, some vague pixels revealed on the left of Fig. 7(a), the second BIMF (see Fig. 7(e)) and residue (see Fig. 7(i)) represent obvious boundary effects. Indeed, it turns out that the other three BEMD methods establish improved results by virtue of various boundary effects alleviation strategies. Nevertheless, there remain ambiguous results in Fig. 7(f) obtained by S-BEMD. Moreover, undesired black blocks appear on the shoulder of the people in
where 1 rx, y rP, ci ðx,yÞ is the ith BIMF, I is supposed to be the total number of BIMFs, and the residue r I ðx,yÞ is taken as an additional element followed by [1]. Consequently, the square of the image IM original ðx,yÞ can be formed as !2 I X 2 ci ðx,yÞ IM original ðx,yÞ ¼ i¼1
¼
I X i¼1
c2i ðx,yÞ þ 2
I I X X
ci ðx,yÞcj ðx,yÞ
ð16Þ
i ¼ 1 j ¼ 1,jai
In case the decomposition result is ideal, the BIMFs will be orthogonal to each other and the cross terms from last part of Eq. (16) should be zero. Regarding these facts, we define the bi-dimensional index of orthogonality (BIO) as below: ! I I P X P X X X 9ci ðx,yÞcj ðx,yÞ9 BIO ¼ ð17Þ 2 i ¼ 1 j ¼ 1,jai x ¼ 1 y ¼ 1 IM original ðx,yÞ where 9 9 denotes the absolute, ci ðx,yÞcj ðx,yÞ gives the dot product between ci ðx,yÞ and cj ðx,yÞ, and IM2original ðx,yÞ presents the dot product of IM original ðx,yÞ and itself. In particular, the BIO reflects the overall orthogonality between BIMFs, which enables us to implement quantitative comparison of the above-mentioned BEMD methods. Moreover, it is worth noting that the nearer the BIO approximates to zero, the better the orthogonality of BIMFs. Table 2 demonstrates the BIO obtained from different BEMD methods. Since there are no boundary effects
132
Z. He et al. / Signal Processing 93 (2013) 124–138
Fig. 6. Boundary extension results of various BEMD methods for 256 256 standard images. The first column: original images, from top to bottom are ‘‘LENA’’, ‘‘BARBARA’’, ‘‘RADISH’’, ‘‘PEPPERS’’, ‘‘TILE’’ and ‘‘FLOOR’’. The second to fourth column: extended results of S-BEMD, T-BEMD and D-BEMD, respectively.
mitigation strategies in O-BEMD, evident boundary effects appear in the BIMFs and the border pixels tend to infinity, which leads to infinite BIO. That means, the BIO of
O-BEMD are invalid. Consequently, we compare the other three BEMD methods in aspect of the orthogonality. By analyzing Table 2, we have to remark that, the BIO
Z. He et al. / Signal Processing 93 (2013) 124–138
133
Fig. 7. Decomposition results of different BEMD methods for ‘‘LENA’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
Fig. 8. Decomposition results of different BEMD methods for ‘‘BARBARA’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
134
Z. He et al. / Signal Processing 93 (2013) 124–138
Fig. 9. Decomposition results of different BEMD methods for ‘‘RADISH’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
Fig. 10. Decomposition results of different BEMD methods for ‘‘PEPPERS’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
Z. He et al. / Signal Processing 93 (2013) 124–138
135
Fig. 11. Decomposition results of different BEMD methods for ‘‘TILE’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
Fig. 12. Decomposition results of different BEMD methods for ‘‘FLOOR’’. The first column: O-BEMD. The second column: S-BEMD. The third column: T-BEMD. The fourth column: D-BEMD.
136
Z. He et al. / Signal Processing 93 (2013) 124–138
Appendix A. Proof of Theorem 1
Table 2 BIO of different BEMD methods. Standard images
LENA BARBARA RADISH PEPPERS TILE FLOOR
Different BEMD methods O-BEMD
S-EMD
T-BEMD
D-BEMD
1 1 1 1 1 1
0.0304 0.0116 0.0681 0.0361 0.0714 0.0576
0.0219 0.0159 0.0530 0.0227 0.0903 0.0532
0.0194 0.0110 0.0481 0.0211 0.0693 0.0455
In this appendix we indicate the proof of Theorem 1. To prove Theorem 1, we need first to discretize the PN1 ð1Þ ð1Þ C-GMC(1,N). Let cðtÞ ¼ dx1 ðtÞ=dt ¼ axð1Þ i ¼ 1 bi xi þ 1 1 ðtÞ þ ðtÞ þ bN , we derive Z k Z k ð1Þ dx1 ðtÞ cðtÞ dt ðA:1Þ dt ¼ dt k1 k1 which can be simplified to Z k ð1Þ ðkÞx ðk1Þ ¼ cðtÞ dt xð1Þ 1 1
ðA:2Þ
k1
obtained from T-BEMD is less than S-BEMD except the images are symmetric (e.g. ‘‘TILE’’) or complicated (e.g. ‘‘BARBARA’’), while D-BEMD deserves the smallest BIO in the three BEMD methods. As discussed above, the BIMFs of D-BEMD involves best orthogonality, while T-BEMD is next, and S-BEMD is worst. In a nutshell, the proposed D-BEMD achieves better performance than the other existing methods, including O-BEMD, S-BEMD, as well as T-BEMD.
Employing the trapezoid formula to Eq. (A.2) can be written explicitly as
5. Conclusions
xð1Þ 1 ðkÞ ¼
Aiming at mitigating boundary effects of BEMD, we have developed a novel boundary extension algorithm inspired by the discrete multivariate gray model. By theoretically proving the equivalence between the continuous and discrete multivariate gray model in Theorem 1, the D-GMC(1,3) is adopted to predict the inner pixels of each extended block in terms of existing border, whose relative data series and characteristic data series are obtained from the coordinates and pixels of the image, respectively. In comparison with existing techniques, including O-BEMD, S-BEMD, as well as T-BEMD, the contribution of our proposed D-BEMD can be outlined in two aspects. On one hand, the boundaries extended by Algorithm 1 accurately reflects the developing trend of original images and yields significant improvements compared with symmetrical extension or texture synthesis based techniques. On the other hand, as stated in Sections 4.2 and 4.3, D-BEMD establishes qualitative and quantitative improvements within reasonable time compared to other three BEMD methods mentioned above. Future research can be conducted, for instance, on the refinement of D-GMC(1,3) to improve the effectiveness of the D-BEMD, the establishment of the algorithmic evaluation system, the decomposition and de-noising of noise images, etc.
Acknowledgment The authors gratefully acknowledge the financial support of the National Basic Research Program of China (973 Program) (2012CB720000), the National Natural Science Foundation of China (Nos. 60975009 and 61171197), IBRSEM, HIT.KLOF.2010017 and Research Fund for the Doctoral Program of Higher Education of China (Nos. 20092302110037 and 20102302110033).
Rk
k1
cðtÞ dt,
1 x1ð1Þ ðkÞx1ð1Þ ðk1Þ ¼ ½cðkÞ þ cðk1Þ 2 " # 1 N 1 X 1 NX ð1Þ ð1Þ ð1Þ ¼ bi xið1Þ ðkÞ þ b x ðk1Þ þ 2b ax ðkÞax ðk1Þ N i iþ1 1 1 þ1 2 i¼1 i¼1
ðA:3Þ Consequently, we can deduce from Eq. (A.3) that N 1 X 10:5a ð1Þ bi x1 ðk1Þ þ xð1Þ i þ 1 ðk1Þ 1þ 0:5a 2ð1 þ 0:5aÞ i¼1
þ xð1Þ i þ 1 ðkÞÞþ
bN 1þ 0:5a
ðA:4Þ
It is interesting to see that Eq. (6) and Eq. (A.4) appear identical form, at this point, we assume the following expression: 8 10:5a > > < a1 ¼ 1 þ 0:5a ðA:5Þ bi > > : bi ¼ , i ¼ 1,2, . . . ,N 1þ 0:5a As previously mentioned, the predicted series of C-GMC(1,N) is evaluated as bN ð1eaðk1Þ Þ a Z k N 1 X þ eaðktÞ bi xð1Þ ðtÞ dt iþ1
ð1Þ aðk1Þ x^ 1 ðkÞ ¼ xð1Þ þ 1 ð1Þe
1
i¼1
aðk1Þ ¼ xð1Þ þ 1 ð1Þe
þ
N1 X i¼1
bi
Z
k 1
bN ð1eaðk1Þ Þ a
eaðktÞ xð1Þ i þ 1 ðtÞ dt,
k ¼ 2,3, . . . ,n þm ðA:6Þ
As referred to [21], the background value of the associated series xð1Þ ðtÞ, t 2 ½k1,k, k ¼ 2,3, . . . ,n þ m, iþ1 i ¼ 1,2, . . . ,N1 can be taken as 12 ðxð1Þ ðk1Þ þxð1Þ ðkÞÞ, iþ1 iþ1 According to Eq. (A.6), it turns out that ð1Þ aðk1Þ x^ 1 ðkÞ xð1Þ þ 1 ð1Þe
þ
N 1 X i¼1
bi
k2 X 1
2 j¼0
bN ð1eaðk1Þ Þ a
ðkjÞ þ xð1Þ ðkj1ÞÞ ðxð1Þ iþ1 iþ1
aðk1Þ ¼ xð1Þ þ 1 ð1Þe
Z
kj
eaðktÞ dt
kj1
N 1 k2 X X bN ð1ea Þeaj ð1eaðk1Þ Þ þ bi a 2a i¼1 j¼0
Z. He et al. / Signal Processing 93 (2013) 124–138
137
Table A1 The relationship between a and z. Values of a Values of z
0.1 0.000075
0.2 0.000533
0.3 0.001575
0.4 0.003200
0.5 0.005208
Values of a Values of z
0.6 0.007200
0.7 0.008575
0.8 0.008533
0.9 0.006075
1 0
ð1Þ ðxð1Þ i þ 1 ðkjÞ þxi þ 1 ðkj1ÞÞ,
k ¼ 2,3, . . . ,n þ m
ðA:7Þ
On the other hand, the predicted series of D-GMC(1,N) can be deduced from ð1Þ x^ 1 ðkÞ
¼ a
ð1Þ 1 x1 ðk1Þ þ
N 1 X
bi 2
i¼1
¼ a1 a1 xð1Þ 1 ðk2Þ þ
i¼1
þ
N1 X
bi
i¼1
2
bi 2
¼ a1 a1 ða
ð1Þ 1 x1 ðk3Þ þ
! ðk2Þ þxð1Þ ðk1ÞÞ þ bN ðxð1Þ iþ1 iþ1
N1 X i¼1
þ
N1 X
bi
i¼1
þ
N1 X
2
bi
i¼1
2
bi 2
ðk3Þ þxð1Þ ðk2ÞÞ þ bN Þ ðxð1Þ iþ1 iþ1
þ
þ
þ
ða1 Þj
ða1 Þj bN ,
bi 2
ðkjÞ þxð1Þ ðkj1ÞÞ ðxð1Þ iþ1 iþ1
k ¼ 2,3, . . . ,nþ m
ðA:8Þ
j¼0
The Maclaurin expansion of ea from Eq. (A.7) and a1 from Eq. (A.8), respectively, yields ea ¼ 1a þ
a2 a3 an þ þð1Þn þ oðan Þ 2! 3! n!
ðA:9Þ
10:5a 1þ 0:5a a a2 an ¼ 1a 1 þ 2 þ þ ð1Þn n þ oðan Þ 2 2 2
a1 ¼
¼ 1a þ
a2 a3 an þ 1 þ þ ð1Þn þ 1 n þ oðan þ 1 Þ 2 22 2
ðA:10Þ
With regard to the small value of a, we calculate the difference between a1 and ea approximately as
z ¼ ða1 Þea
a3 23
þ
a3 a4 a4 a3 þa4 ¼ þ 3! 23 4! 12
ðA:11Þ
from which we obtain the numerical relationship between a and z (see Table A1). As a consequence, there exists a1 ea , which leads to bi b 10:5a b b ¼ i ½1ða1 Þ i 1ea , ¼ i 1 bi ¼ 1 þ0:5a 1 þ 0:5a a a a i ¼ 1,2, . . . ,N
ðA:12Þ
Note that a1 e , bi ðbi =aÞð1e Þ, i ¼ 1,2, . . . ,N Pk2 aj and ð1ea Þ ¼ 1eaðk1Þ , Eq. (A.8) can be j¼0e a
ð1Þ ðxð1Þ i þ 1 ðkjÞ þ xi þ 1 ðkj1ÞÞ
k2 X bN
a
bi ð1Þ ð1ea Þeaj ðxð1Þ i þ 1 ðkjÞ þ xi þ 1 ðkj1ÞÞ 2a
eaj ð1ea Þ
aðk1Þ ¼ xð1Þ 1 ð1Þe
ðk1Þ þxð1Þ ðkÞÞ þ bN ðxð1Þ iþ1 iþ1
a
N 1 X k2 X
þ
bi ð1ea Þeaj ðxð1Þ ðkjÞ þ xð1Þ ðkj1ÞÞ iþ1 iþ1 2a i¼1j¼0
þ
bN ð1eaðk1Þ Þ a
¼ ða
k2 X
N 1 X k2 X
j¼0
N 1 X k2 X
2
ða1 Þj bN
i¼1j¼0
k1 ð1Þ x1 ð1Þ 1Þ
þ
k2 X
bi
j¼0
# ðk2Þ þxð1Þ ðk1ÞÞ þ bN ðxð1Þ iþ1 iþ1
i¼1j¼0
ða1 Þj
aðk1Þ xð1Þ 1 ð1Þe
¼
þ
N 1 X k2 X i¼1j¼0
ðk1Þ þxð1Þ ðkÞÞ þ bN ðxð1Þ iþ1 iþ1
"
ð1Þ x^ 1 ðkÞ ¼ ða1 Þk1 xð1Þ 1 ð1Þ
þ
ðk1Þ þxð1Þ ðkÞÞ þ bN ðxð1Þ iþ1 iþ1 N1 X
rewritten as
ðA:13Þ
Up till now, it is notable that Eq. (A.7) is the same as Eq. (A.13), which verifies the equivalence between C-GMC(1,N) and D-GMC(1,N). This completes the outline of proof of Theorem 1. References [1] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, E.H. Shih, Q. Zheng, C.C. Tung, H.H. Liu, The empirical mode decomposition method and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proceedings of the Royal Society of London, Series A (1998) 903–995. [2] Z.G. Xu, B.X. Huang, F. Zhang, Improvement of empirical mode decomposition under low sampling rate, Signal Processing 89 (11) (2009) 2296–2303. [3] N. Tsakalozos, K. Drakakis, S. Rickard, A formal study of the nonlinearity and consistency of the empirical mode decomposition, Signal Processing 92 (9) (2012) 1961–1969. [4] J.-C. Nunes, Y. Bouaoune, E. Delechelle, O. Niang, P. Bunel, Image analysis by bidimensional empirical mode decomposition, Image and Vision Computing 21 (12) (2003) 1019–1026. [5] S.M.A. Bhuiyan, N.O. Attoh-Okine, K.E. Barner, A.Y. Ayenu-Prah, R.R. Adhami, Bidimensional empirical mode decomposition using various interpolation techniques, Advances in Adaptive Data Analysis 1 (2) (2009) 309–338. [6] Z.H. Wu, N.E. Huang, X.Y. Chen, The multi-dimensional ensemble empirical mode decomposition method, Advances in Adaptive Data Analysis 1 (3) (2009) 339–372. [7] M.-H. Yeh, The complex bidimensional empirical mode decomposition, Signal Processing 92 (2) (2012) 523–541. [8] J. Fleureau, A. Kachenoura, L. Albera, J.-C. Nunes, L. Senhadji, Multivariate empirical mode decomposition and application to multichannel filtering, Signal Processing 91 (12) (2011) 2783–2792. [9] Y. Xu, B. Liu, J. Liu, S. Riemenschneider, Two-dimensional empirical mode decomposition by finite elements, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science (2006) 3081–3096.
138
Z. He et al. / Signal Processing 93 (2013) 124–138
[10] S.M.A. Bhuiyan, R.R. Adhami, J.F. Khan, Fast and adaptive bidimensional empirical mode decomposition using order-statistics filter based envelope estimation, EURASIP Journal on Advances in Signal Processing 2008 (21) (2008) 1–18. [11] R. Pabel, R. Koch, G. Jager, A. Kunoth, Fast empirical mode decompositions of multivariate data based on adaptive spline-wavelets and a generalization of the Hilbert–Huang-transform (HHT) to arbitrary space dimensions, Advances in Adaptive Data Analysis 2 (3) (2010) 337–358. [12] A. Linderhed, Adaptive image compression with wavelet packets and empirical mode decomposition, Department of Electrical Engi¨ neering, Linkoping University, 2004. [13] A. Linderhed, Image empirical mode decomposition: a new tool for image processing, Advances in Adaptive Data Analysis 1 (2) (2009) 265–294. [14] Y. Tian, K. Zhao, Y.P. Xu, F.Y. Peng, An image compression method based on the multi-resolution characteristics of BEMD, Computers & Mathematics with Applications 61 (8) (2011) 2142–2147. [15] P.-H. Chen, Y.-C. Yang, L.-M. Chang, Illumination adjustment for bridge coating images using BEMD-Morphology approach (BMA), Automation in Construction 19 (4) (2010) 475–484. [16] J.N. Huang, B.B. Zhao, Y.Q. Chen, P.D. Zhao, Bidimensional empirical mode decomposition (BEMD) for extraction of gravity anomalies associated with gold mineralization in the Tongshi gold field,
[17]
[18]
[19]
[20]
[21] [22] [23]
[24]
Western Shandong Uplifted Block, Eastern China, Computers & Geosciences 36 (7) (2010) 987–995. S.M.A. Bhuiyan, J.F. Khan, R.R. Adhami, A novel approach of edge detection via a fast and adaptive bidimensional empirical mode decomposition method, Advances in Adaptive Data Analysis 2 (2) (2010) 171–192. Z.X. Liu, S.L. Peng, Boundary processing of bidimensional EMD using texture synthesis, IEEE Signal Processing Letters 12 (1) (2005) 33–36. A.A. Efros, T.K. Leung, Texture synthesis by non-parametric sampling, in: Proceedings of the Seventh IEEE International Conference on Computer Vision, 1999, pp. 1033–1038. N.-M. Xie, S.-F. Liu, Discrete grey forecasting model and its optimization, Applied Mathematical Modelling 33 (2) (2009) 1173–1186. T.-L. Tien, A research on the grey prediction model GM(1,N), Applied Mathematics and Computation 218 (9) (2012) 4903–4916. T.-L. Tien, The indirect measurement of tensile strength by the new model FGMC (1,N), Measurement 44 (10) (2011) 1884–1897. Z. He, Y. Shen, Q. Wang, Boundary extension for Hilbert–Huang transform inspired by gray prediction model, Signal Processing 92 (3) (2012) 685–697. ¨ V. Kwatra, A. Schodl, I. Essa, G. Turk, A. Bobick, Graphcut textures: image and video synthesis using graph cuts, ACM Transactions on Graphics 22 (3) (2003) 277–286.