A green’s function-based Bi-dimensional empirical mode decomposition

Highlights

• A new GiT-BEMD algorithm is developed.

• Based on Green’s function with tension.

• Much faster than canonical BEMD algorithms.

• Easy to implement and very robust.

• Superior performance tested on artificial and real images.

Abstract

Bidimensional Empirical Mode Decomposition(BEMD) interprets an image as a superposition of Bidimensional Intrinsic Mode Functions (BIMFs). They are extracted by a process called sifting, which encompasses two-dimensional surface interpolations connecting a set of local maxima or minima to form corresponding envelope surfaces. Existing surface interpolation schemes are computationally very demanding and often induce artifacts in the extracted modes. This paper suggests a novel method of envelope surface interpolation based on Green’s functions. Including surface tension greatly improves the stability of the new method which we call Green’s function in tension-based BEMD (GiT-BEMD). Simulation results, using toy images with various textures, facial images and functional neuroimages, demonstrate the superior performance of the new method when compared to its canonical BEMD counterpart. GiT-BEMD strongly speeds up computations and achieves a higher quality of the extracted BIMFs. Furthermore, GiT-BEMD can be extended simply to an ensemble-based variant (GiT-BEEMD), if needed. In summary, the study suggests the new variant GiT-BEMD as a highly competitive, fast and stable alternative to existing BEMD techniques for image analysis.

Introduction

Empirical mode decomposition (EMD), as pioneered by Huang et al. [21], is a data driven signal processing algorithm that quantitatively decomposes any nonlinear and non-stationary data into intrinsic modes, thereby obtaining local features, and their related time-frequency distribution. The first step of this method decomposes the data/signal into its characteristic intrinsic mode functions (IMFs) [12], [58] while the second step finds the time frequency distribution of the data from each IMF by utilizing the concepts of Hilbert transform and instantaneous frequency. The complete process is also known as the Hilbert–Huang Transform (HHT) [4]. Lacking any rigorous mathematical basis, [68] advocated local physical rather than global mathematical constraints to preserve any physical meaning of the IMFs extracted. Soon after its invention, this decomposition technique has been extended to analyze multi-dimensional data sets [17], [36], [37], [42], [67] including complex-valued data sets [4], [23], [33], [56] and implementations on GPUs for parallel processing [11]. Besides an extension to multi-dimensional data sets, EMD also has been extended to multi-variate data sets, most notably multi-channel recordings of biomedical signals [18], [24], [40], [43], [44], [45], [46], [50]. Since its invention numerous applications of EMD have been reported in such diverse fields as functional neuroimaging [19], [41], [49], face recognition [20], facial emotion recognition [3], biomedical signals [5], [25], [35], [39], [47], [52], [71], [72], neuromonitoring [16], analysis of complex networks [27], image enhancement [10], fault diagnosis of mechanical systems [26], ultrasound echo detection [34], speaker identification [65], speech enhancement [24], forecasting [1], moving target recognition [73] etc. just to mention a few more recent publications. However, most of these applications concern one-dimensional data sets and corresponding EMD algorithms.

Obviously, two-dimensional image data sets were of special interest [48]. In a first approach, such two-dimensional data was treated as a collection of one-dimensional slices, which were decomposed with one-dimensional EMD [31], [32]. This procedure is called pseudo-two-dimensional EMD [67]. The latter technique treats each row and/or each column of the 2D data set separately by a 1D EMD, which renders the sifting process faster than in a genuine 2D decomposition. But this parallel 1D implementation results in poor BIMF components compared to the canonical 2D procedure due to the fact that the former ignores the correlation among the rows and/or columns of a 2D image [29]. In addition, pseudo-2D-EMD needs a coherence structure associated with the spatial scales in a particular direction, which significantly limits its use. These recent developments in analysis methods for non-linear and non-stationary data sets have received considerable attention by image analysts. Thus several attempts have been started lately to extend EMD to multi-dimensional data sets like two-dimensional (2D) data arrays and images. These extensions are variously known as bidimensional EMD (BEMD), image EMD (IEMD), 2D EMD and so on [13], [28], [30], [31], [36], [37], [38], [69], [70]. Some of these works especially exploit mode decompositions to compute texture information contained in the images. In [74] a new two-dimensional EMD (2DEMD) method is proposed, which is claimed being faster and better-performing than current 2DEMD methods. In [69] rotation-invariant texture feature vectors are extracted at multiple scales or spatial frequencies based on a BEMD. In Nunes et al. [36], [37], [38] the BEMD-based texture extraction algorithm is demonstrated in experiments with both artificial and natural images. A major breakthrough has been achieved by Wu et al. [67], who recently proposed a Multi-dimensional Ensemble Empirical Mode Decomposition (mdEEMD) for multidimensional data arrays. This algorithm turned out to be very efficient in practical two-dimensional applications, especially if combined with Ensemble Empirical Mode Decomposition (EEMD) [66]. Lately also a full Bayesian approach, based on a reversible jump Markov Chain Monte Carlo procedure to sample from the unknown posteriors, has been proposed by Bouguila and Elguebaly [8] and applied for image texture retrieval and classification. Although this represents an optimal data analysis technique, its computational complexity is prohibitive, especially when big data have to be considered.

The computationally most demanding operation of all these algorithms involves an envelope surface interpolation step. While cubic spline interpolation is preferred for 1D interpolation [57], various types of radial basis function, multilevel B-spline, Delaunay triangulation, Order-Statistics Filter, Finite Elements method and so on have been used for 2D scattered data interpolation [21], [70]. Among them, Delaunay triangulation and a Finite Elements Method provide relatively faster decomposition compared to the other methods. In [13], the influence of various interpolation methods is studied, and a sifting process is proposed based on a Delaunay triangulation with subsequent cubic interpolation on triangles. Subsequently, the envelope surface interpolation step is replaced by either a direct envelope surface estimation method or radial basis function interpolators [6], [7]. Finally, the modified 2D EMD algorithm proposed in [69] implements the FastRBF algorithm for the estimation of the envelope surfaces.

However, in multi-dimensional EMD it is generally required finding local maxima and local minima, jointly known as local extrema, and subsequently interpolating these points in each iteration of the process. In general, surface interpolation schemes have to face the following issues:

Computational load: Local extrema of an one-dimensional (1D) signal are obtained using either a sliding window or local derivative, and local extrema of a 2D data/image set are extracted using either a sliding window technique or various morphological operations [21], [36], [37], [38]. Hence, detection and interpolation of local extrema during each iteration render the process complicated and time consuming. The situation is more difficult for the case of BEMD as it requires an interpolation of a set of 2D scattered two-dimensional data arrays during each iteration. For some images decomposition may take hours or days unless any additional stopping criterion is employed. However, additional stopping criteria may result in an inaccurate and incomplete decomposition [22], [53].

Boundary artifacts: Another common and significant problem related to the interpolation of scattered 2D data in BEMD is that the maxima or minima map often does not contain any data points (interpolation centers) at the boundary region. This shortcoming becomes especially severe for bidimensional intrinsic modes, henceforth called BIMFs, extracted at the end of the decomposition process where spatial frequencies are low and the modes encompass only few extremal points. Currently available interpolation methods for scattered data are inefficient in handling this kind of situation. Additionally, the effect of incorrect interpolation at the boundary gradually propagates towards the central region from iteration to iteration, and from BIMF mode to BIMF mode, causing corrupted BIMFs.

Over-and Undershooting: Overshooting or undershooting is another problem of interpolation-based envelope estimation, which causes incorrect BIMFs. This is especially true for higher order spline interpolation schemes.

Although a few modifications have been suggested in the literature to reduce the number of iterations and/or to overcome boundary artifacts [13], [22], [29], the technique still suffers from the above-mentioned problems to some extent. In the BEMD process, the number of extrema decreases from BIMF to BIMF. For the later modes, there may be very few irregularly spaced local maxima or minima, which can cause highly erroneous and misleading upper or lower envelopes, and thus result in highly incorrect BIMFs. In order to improve the algorithm performance, some modifications have been suggested for one-dimensional EMD [54], [59], which may not be useful for BEMD in the context of processing speed and algorithm complexity. For example, in fast and adaptive BEMD (FABEMD) [6], although it enhances the computational load noticeably with artificial texture, sometimes it may impose some difficulty in applying the FABEMD algorithm with real images; because it is difficult to determine the proper type of this method without a priori knowledge, which leads to improper intrinsic modes. On the other hand, the requirement for manipulation of window width for a BIMF may impose additional complexity, when the calculated value does not appear larger than the previous BIMF mode. Due to the property of order statistics, filter-based envelope estimation followed by a smoothing operation, oversifting in FABEMD may cause improper BIMFs as well. Moreover, any type of additional processing steps may turn the process more complex and computationally expensive.

Nonetheless, BEMD is a promising image processing technique that has been applied successfully in various real world problems, for example, medical image analysis [2], image analysis [36], texture analysis [30], and so on. But the main problem for BEMD is either a high computational load or otherwise a reduced quality of the extracted intrinsic modes. Hence, any improvement of any BEMD algorithm is very important. Following we refer to the BEEMD approach developed by Wu et al. [67] as the canonical BEMD variant, and contrast it with our novel bi-dimensional EMD (GiT-BEMD) approach that replaces the direct spline interpolation step by an interpolation with splines in tension employing Green’s functions [62]. Although, cubic splines are in widespread use, because of their smooth shape, these functions can exhibit unwanted oscillations between data points. Adding tension to the spline overcomes this drawback. Here, we apply a technique for interpolation and gridding in bi-dimensional EMD applications using Green’s functions for splines in tension, and examine some of the properties of these functions. The technique is borrowed from geophysics where it is in use already since more than a decade [61]. Physical sciences have a frequent need for data interpolation and gridding. Such tasks are commonly accomplished by averaging [60] and finite difference methods [9] employing cubic splines on a regular grid. Introducing surface tension often helps to suppress undesired oscillations of such splines [55]. The Generic Mapping Tools offer a software package implementing an algorithm which uses continuous curvature splines in tension [63], [64]. Minimum curvature gridding based on Green’s functions of the bi-harmonic operator has been proposed first by Sandwell [51]. It offered enhanced flexibility by employing both data values and gradients to constrain the interpolating surface. The method further allowed least squares fitting to noisy data sets and could evaluate the surface at any location instead of being confined to a regular grid. However, the appearance of extraneous inflection points, common to all minimum curvature methods, still represented a major obstacle to applications. Wessel and Bercovici [62] generalized the approach of Sandwell by including surface tension to the Green’s function expansion of the interpolating surface. In summary, for moderate amounts of data, the Green’s function technique is superior to conventional finite difference methods because both data values and directional gradients can be used to constrain the model surface. Also, noise can be suppressed easily by striving a least-squares fit rather than considering a strict interpolation, and the model can be estimated at arbitrary locations rather than only on a rectangular grid [62]. Moreover, including surface tension greatly improves the stability of the method relative to gridding without tension. Recently, Wessel presented a Green’s function based general purpose interpolator for both Cartesian and spherical surface data called Greenspline [61].

In this study we propose, for the first time, to join the Green’s function method with the BEMD technique thus eliminating the poor interpolation effects and reducing the computation time for each iteration. Most favorably, this interpolation technique only needs very few iterations for estimating each BIMF as it does not need to employ an ensemble. The proposed fast and stable GiT-BEMD method thus can be a good alternative providing an efficient BEMD processing.

The organization of the paper is as follows: before introducing the novel concept of GiT-BEMD, the regular BEMD process is briefly summarized in Section 2 of this paper. The detailed description of the proposed GiT-BEMD algorithm is given in Section 3. Although the method of detecting local extrema suggested in GiT-BEMD is the same as in normal BEMD, it is explained in the first part of Section 3 to further understanding of the proposed envelope estimation technique, since it requires information about local extrema as its basic ingredients. The second part of Section 3 describes the new method of envelope estimation. Simulation results with various images comparing GiT-BEMD and canonical BEEMD are given in Section 4. Finally, a conclusion is drawn in Section 5.