Copyright © 2008 Sharif M. A. Bhuiyan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

A novel approach for bidimensional empirical mode decomposition (BEMD) is proposed in this paper. BEMD decomposes an image into multiple hierarchical components known as bidimensional intrinsic mode functions (BIMFs). In each iteration of the process, two-dimensional (2D) interpolation is applied to a set of local maxima (minima) points to form the upper (lower) envelope. But, 2D scattered data interpolation methods cause huge computation time and other artifacts in the decomposition. This paper suggests a simple, but effective, method of envelope estimation that replaces the surface interpolation. In this method, order statistics filters are used to get the upper and lower envelopes, where filter size is derived from the data. Based on the properties of the proposed approach, it is considered as fast and adaptive BEMD (FABEMD). Simulation results demonstrate that FABEMD is not only faster and adaptive, but also outperforms the original BEMD in terms of the quality of the BIMFs.

1. Introduction

Empirical mode decomposition (EMD), originally developed by Huang et al. [1, 2], is a data driven signal processing algorithm that has been established to be able to perfectly analyze nonlinear and nonstationary data by obtaining local features and time-frequency distribution of the data. The first step of this method decomposes the data/signal into its characteristic intrinsic mode functions (IMFs), 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) [1]. This decomposition technique has also been extended to analyze two-dimensional (2D) data/images, which is known as bidimensional EMD (BEMD), image EMD (IEMD), 2D EMD and so on [3–8]. Both EMD and BEMD require finding local maxima and local minima points (jointly known as local extrema points) and subsequent interpolation of those points in each iteration of the process. Local extrema points of one-dimensional (1D) signal are obtained using either a sliding window or local derivative, and local extrema points of 2D data/image are extracted using sliding window or various morphological operations [1–4]. Cubic spline interpolation is preferred for 1D interpolation while various types of radial basis function, multilevel B-spline, Delaunay triangulation, finite-element method, and so on have been used for 2D scattered data interpolation [1–7], where Delaunay triangulation and finite-element method provide relatively faster decomposition compared to the other methods. Beside 2D implementation of the BEMD process, 1D EMD has also been applied to images to obtain 2D IMFs or bidimensional IMFs (BIMFs) [8–10]. In this technique, each row and/or each column of the 2D data is processed by 1D EMD, which makes it a faster process. However, it has been found that this 1D implementation results in poorer BIMF components compared to the standard 2D procedure due to the fact that the former ignores the correlation among the rows and/or columns of a 2D image [11].

In EMD or BEMD, extraction of each IMF or BIMF requires several iterations. Hence, extrema detection and interpolation at each iteration make the process complicated and time consuming. The situation is more difficult for the case of BEMD that requires 2D scattered data interpolation at each iteration. For some images it may take hours or days for decomposition unless any additional stopping criterion is employed, whereas additional stopping criteria may result in inaccurate and incomplete decomposition [12–15] that may not be desired. Another common and significant problem related to the 2D scattered data interpolation in BEMD is that the maxima or minima map often does not contain any data points (interpolation centers) at the boundary region, which may be more severe for the later modes of decomposition. Currently available scattered data interpolation methods are inefficient in handling this kind of situation. Additionally, the effect of incorrect interpolation at the boundary gradually propagates into the mid region from iteration to iteration and from BIMF mode to BIMF mode causing corrupted BIMFs. Overshooting or undershooting is another problem of interpolation-based envelope estimation, which causes incorrect BIMFs. Although a few modifications have been suggested in the literature to reduce the number of iterations and/or to overcome the boundary effects [6, 11, 12], the technique still suffers from the above-mentioned problems to some extent. In the BEMD process, the number of extrema points decreases from one mode to the next mode. For the later modes, there may be very few irregularly spaced local maxima or minima points, which can cause highly erroneous and misleading upper or lower envelopes, and thus incorrect modes of BIMFs. In order to improve the algorithm performance, some modifications have been suggested for EMD [16–20], which may not be useful for BEMD in the context of processing speed and algorithm complexity. Moreover, any types of additional processing steps may make the process more complex and computationally extensive.

BEMD is a promising image processing algorithm that can be applied successfully in various real world problems, for example, medical image analysis, pattern analysis, texture analysis, and so on. But the problem, due to scattered data interpolation in BEMD, limits its application to very small size images, while the size of the real images may be much bigger than is suitable for BEMD processing. It is also not appropriate to reduce the size of the images only for the purpose of BEMD processing and thus loose the fine details and/or relevant information. Hence, improvement of the BEMD algorithm is very important. In this paper, a novel BEMD approach is suggested that replaces the interpolation step by a direct envelope estimation method. In this technique, spatial domain sliding order-statistics filters, namely, MAX and MIN filters, are employed to get the running maxima and running minima of the data, which is followed by smoothing operation to get the upper envelope and lower envelope, respectively. The size of the order-statistics filters is derived from the available information of maxima and minima maps. In addition to eliminating the poor interpolation effects and reducing the computation time for each iteration, this process facilitates performing only one iteration for each BIMF. The proposed fast and adaptive BEMD (FABEMD) method can be a good alternative for efficient BEMD processing.

For ease of discussion, some new terms have been introduced in this paper in place of the existing terms associated with EMD or BEMD. Before introducing the novel concepts of FABEMD, the regular BEMD process is briefly reviewed in Section 2 of this paper. The detailed description of the proposed FABEMD algorithm is given in Section 3. Although the extrema detection method suggested in FABEMD is the same as in BEMD, it is explained in the first part of Section 3 for understanding the proposed envelope estimation technique, since it requires the extrema information as its foundation. The second part of Section 3 describes the new method of envelope estimation. Simulation results with various images comparing FABEMD and BEMD are given in Section 4. Finally, concluding remarks are given in Section 5.

2. BEMD Overview

EMD or BEMD is a sifting process that decomposes a signal into its IMFs or BIMFs and a residue based basically on the local frequency or oscillation information. The first IMF/BIMF contains the highest local frequencies of oscillation or the highest local spatial scales, the final IMF/BIMF contains the lowest local frequencies of oscillation and the residue contains the trend of the signal/data. Like time-frequency distribution with EMD, acquiring the space-spatial-frequency distribution of 2D data/image is possible with BEMD, which may be named as bidimensional HHT (BHHT). Although direct estimation of the horizontal and vertical frequencies of BIMFs has been studied [21], BHHT has not yet been reported in the literature. It is claimed and experimentally shown that the HHT performs better than the other existing techniques of analyzing the time-frequency distribution of nonstationary and nonlinear data [1]. Thus, HHT or BHHT can better represent the local frequency and amplitude scale of the signal if the IMF or BIMF components appear perfect. However, decomposition of an image into BIMFs alone can offer a wide variety of image processing applications. Hence, the following discussion will be limited to the first part of BEMD only, that is, decomposition of an image into BIMFs and the Residue. It should be noted that once the BIMFs are obtained, the space-spatial-frequency distribution of an image can be acquired with standard techniques of 2D Hilbert spectral analysis (HSA).

2.1. Properties of IMF/BIMF

The IMFs of a signal obtained by EMD are expected to have the following properties [1, 2, 12, 22].

(i) In the whole data set, the number of local extrema (maxima and minima together) and the number of zero crossings must be equal or differ by at most one. (ii) There should be only one mode of oscillation, that is, only one local maxima or local minima, between two successive zero crossings. (iii) At any point, the mean value of the upper and lower envelopes, defined by the local maxima and minima points, is zero or nearly zero. (iv) The IMFs are locally orthogonal among each other and as a set. In fact, property (i) ensures property (ii), and vise versa. The definition and properties of the BIMFs are slightly different from the IMFs. It is sufficient for BIMFs to follow only the final two (iii) and (iv) properties given above [3, 4]. In fact, due to the properties of an image and the BEMD process, it is not possible to satisfy the first two properties (i) and (ii) given above in the case of BIMFs, since the maxima and minima points are defined in a 2D scenario for an image. For the same reason, it is also difficult or impossible to define and/or to achieve any characteristic relationships between the number of maxima points and the number of minima points for BIMFs.

2.2. Steps of BEMD

The required properties of IMFs are achieved via an “empirical” iterative process [1] in EMD. The same algorithm applies for BEMD as well, where extrema detection and interpolation are carried out using 2D versions of the corresponding 1D methods. Let the original image be denoted as , a BIMF as , and the residue as . In the decomposition process ith BIMF is obtained from its source image , where is a residue image obtained as and . It requires one or more iterations to obtain , where the intermediate temporary state of BIMF (ITS-BIMF) in jth iteration can be denoted as . With the definition of the variables, the steps of the BEMD process can be summarized as follows [1–5].

(i) Set . Take I and set . (ii) Set . Set . (iii) Obtain the local maxima map (LMMAX) of , denoted as . (iv) Form the upper envelope (UE) of , denoted as by interpolating the maxima points in . (v) Obtain the local minima map (LMMIN) of , denoted as . (vi) Form the lower envelope (LE) of , denoted as by interpolating the minima points in . (vii) Find the mean/average envelope (ME) as . (viii) Calculate as . (ix) Check whether follows the BIMF properties. These criteria are verified by finding the standard deviation (SD), denoted as D, between and as defined below and comparing it to the desired threshold [1, 3] (1a) where denotes the coordinate of the 2D data, M is the total number of rows and N is the total number of columns in the 2D data. The SD can also be defined as (1b)Although both of the SD measures in (1a) and (1b) provide a global measure of SD, the later one is not dominated by the local fluctuations of the denominator. Normally, a low value of D (e.g., below 0.5 for (1a) and below 0.05 for (1b)) is chosen to ensure nearly zero envelope mean of the BIMF. (x) If meets the criteria given in step (ix), then take . Set first and then . Go to step (xi). If does not meet the stopping criteria, then set , go to step (iii) and continue up to step (x) as before until the criteria are fulfilled. (xi) Determine whether has less than three extrema points, and if so, this is the residue of the image (i.e., ); and the decomposition is complete. Otherwise, go to step (ii) and continue up to step (xi) to obtain the subsequent BIMFs. In the process of extracting the BIMFs, the number of extrema points in should be lower than that in . The BIMFs and the residue of an image together can be named as bidimensional empirical mode components (BEMCs). Except for the truncation error of the digital computer, the summation of all BEMCs returns the original data/image back as given by (1) where is the ith BEMC and K is the total number of BIMFs excluding the residue. An orthogonality index (OI), denoted as O, has been proposed for IMFs in [1], which may be extended for the case of BEMCs as follows: (2) A low value of OI indicates a good decomposition in terms of local orthogonality among the BEMCs. In general, OI values less than or equal to 0.1 are acceptable.

2.3. Issues Related to BEMD

The decomposition of an image into BEMCs is not a unique process. The number of BEMCs and their characteristics depend on the extrema detection method, interpolation technique, and stopping criteria of the iterations for each BIMF. In that sense, there are an infinite number of BEMC sets for each image [12]. As mentioned in Section 1, local extrema (maxima and minima) points of 2D data/image are obtained using 2D sliding window or various morphological operations [1–4] and radial basis function, multilevel B-spline, Delaunay triangulation, finite-element method, and so forth, 2D scattered data interpolation [3–7] have been used for interpolating the extrema points to form the upper and lower envelopes. To stop the iterations for each BIMF, the SD threshold criterion is mostly used to satisfy the zero envelope mean, although there are several additional stopping criteria that may be employed [12–15]. The performance of scattered data interpolation in the BEMD process is highly dependent on the interpolation centers, their orientation, location, numbers, and so on. Hence, local maxima and minima maps play a significant role in creating the upper and lower envelopes. Absence or lack of extrema points at the boundaries of ITS-BIMFs and the presence of very few extrema points in the source images for higher values of i, cause erroneous surface interpolation that results in misleading upper or lower envelopes and hence incorrect BIMFs. Because the surface interpolation method fits a surface in iterative optimization approach utilizing the scattered data arising from the extrema points, it makes the BEMD process an extremely slow one.

3. FABEMD Algorithm Details

With the intention of overcoming the difficulty in implementing BEMD via the application of surface interpolation, a novel approach is devised that eliminates the need for surface interpolation. This new BEMD process, named as fast and adaptive BEMD (FABEMD), differs from the actual BEMD algorithm, basically in the process of estimating the upper and lower envelopes and in limiting the number of iterations per BIMF to one. Hence, the steps of the FABEMD algorithm remain the same as BEMD given in Section 2.2 with maximum required value of j (iteration index for each BIMF) equal to one considered being sufficient. The details of extrema detection and envelope formation of the FABEMD process are discussed in this section.

3.1. Detection of Local Extrema

Detection of local extrema means finding the local maxima and minima points from the given data. The 2D array of local maxima points is called a maxima map (LMMAX) and the 2D array of local minima points is called a minima map (LMMIN). Like BEMD, neighboring window method is employed to find local maxima and local minima points from the jth ITS-BIMF of any source image , where for . In this method, a data point/pixel is considered as a local maximum (minimum), if its value is strictly higher (lower) than all of its neighbors. Let A be an 2D matrix represented by (3) where a_mn is the element of A located in the mth row and nth column. Let the window size for local extrema determination be . Then, (4) where (5)Generally, a window (i.e., ) results in an optimum extrema map for a given 2D data. However, a higher window size may be used in some applications; but this will result in a lower, if not equal, number of local extrema points for a given data matrix. Let us consider the data matrix given in Figure 1(a) for illustration purposes. The maxima map given in Figure 1(b) and minima map given in Figure 1(c) are obtained when a neighboring window is used for every point in the matrix. For finding extrema points at the boundary or corner, the neighboring points within the window that are beyond the image are neglected. As an example, windows centered at a₃₂, a₇₅, and a₂₆ with darker grids are also shown in Figure 1(a). The center element of the first window is a local maximum, the center element of the second window is a local minimum, while the center element of the third window is neither a local maximum nor a local minimum.

Figure 1: (a) A sample data matrix; (b) local maxima map obtained from (a); and (c) local minima map obtained from (a).

3.2. Generating Upper and Lower Envelopes

After obtaining the maxima and minima maps, and , respectively, from a given ITS-BIMF , the next step is to create the continuous upper and lower envelopes, and . In usual BEMD, suitable 2D scattered data interpolation is applied to and to create these envelopes. In this work, a simple but efficient modification has been formulated for the generation of upper and lower envelopes. This approach basically applies two order statistics filters to approximate the envelopes, where a MAX filter is used for upper envelope and a MIN filter is used for lower envelope. Order statistics filters are spatial filters whose response is based on ordering (ranking) the elements contained within the data area encompassed by the filter [23]. The response of the filter at any point is determined by the ranking result. The crucial part of applying the order statistics filters for envelope estimation is to determine an appropriate size for the filter. Based on the desired properties of BIMFs along with the characteristics of and for a given , the method described in Section 3.2.1 is developed for window size determination to extract the corresponding BIMF .

3.2.1. Determining Window Size for Order-Statistics Filters

The window size for order statistics filters is determined based on the maxima and minima maps obtained from a source image , that is, based on and derived from when and . For each local maximum point in , that is, for each nonzero element in , the Euclidean distance to the nearest nonzero element is calculated. The array of distances thus obtained is called adjacent maxima distance array (AMAXDA), denoted as , where the number of elements in AMAXDA is equal to the number of local maxima points in the maxima map . Figure 2(a) shows the maxima map of Figure 1(b) with the maxima points being represented as bright boxes while the other points are represented as dark boxes. Figures 2(b) and 2(c) show two points of interest from the set of maxima points marked with “” and their corresponding nearest neighbors marked with “”.

Figure 2: (a) Maxima map of Figure 1(b) shown with shades where the brighter boxes represent the location of the maxima points, (b) and (c) sample maxima point and its nearest neighbor shown, respectively, with “” and “”.

Similarly, the array of distances obtained from the local minima map is called adjacent minima distance array (AMINDA), denoted as , where the number of elements in AMINDA is equal to the number of local minima points in the minima map . Both and are sorted in descending order for convenient selection of distances from these arrays. Considering square window, the gross window width for order statistics filters can be selected in many different ways using the distance values in and among which four choices are given below (6) where denotes the maximum value of the elements in the array and denotes the minimum value of the elements in the array . is then rounded to the nearest odd integer to get the final window width producing a window of size . The relation of the distances obtained from (6) is . Let the order statistics filter widths (OSFW) obtained via (6) be defined as Type-1, Type-2, Type-3, and Type-4, respectively, where Type-1 and Type-4 may also be denoted as lowest distance OSFW (LD-OSFW) and highest distance OSFW (HD-OSFW), respectively. required for th BIMF generally appears larger than that for the ith BIMF if using Type-3 or Type-4 OSFW; however, for th BIMF sometimes may not appear larger than that for the ith BIMF if using Type-1 or Type-2 OSFW. Therefore, if the calculated for a BIMF mode is not larger than the previous BIMF mode, then additional manipulation may be required to make it larger than the previous mode (e.g., current may be taken as approximately 1.5 times of the previous ). Though it is not necessary, it will ensure the currently existing properties of BIMF hierarchy in the sense that the later BIMF will contain coarser local spatial scales [1, 3]. It will be clear from Section 4 that the choice of from the above four options depends on the application and/or desired BIMF characteristics.

It is preferable to apply the same window size for both MAX and MIN filters as discussed above, though it may be possible to choose different window sizes for them. For example, window size for the MAX filter can be selected based on the distances in AMAXDA, while window size for the MIN filter can be selected based on the distances in AMINDA as follows: (7)(8) Equation (7) can be used for the MAX filter and (8) can be used for the MIN filter. However, there is a practical limitation to this approach. In some situations, there may be only one local maxima (minima) in a source image , which will result in an empty array for and thus will prevent upper (lower) envelope formation and hinder the algorithm before it satisfies the extrema criteria for stopping. On the other hand, employing the same size for MAX and MIN filters for the same BIMF induces extraction of similar spatial scales into that BIMF, while different window sizes for MAX and MIN filters may obstruct this process. It is worthwhile to mention an additional option for the selection of before describing the envelope formation in Section 3.2.2. Based on the image or desired properties of BIMFs, may be chosen arbitrarily as well. In that case, for th BIMF should be chosen higher than the for the ith BIMF; but extraction of BIMFs will be less data driven with an arbitrary selection of . The various possibilities of window sizes for MAX and MIN filters for envelope formation provide different decomposition of an image. It is this feature that makes the proposed approach an adaptive one.

3.2.2. Applying Order Statistics and Smoothing Filters

With the determination of window size for envelope formation, MAX and MIN filters are applied to the corresponding ITS-BIMF to obtain the upper and lower envelopes, and , as specified below: (9)(10)In (9) the value of the upper envelope at any point is simply the maximum value of the elements in in the region defined by , where is the square region of size centered at any point of . Similarly, in (10) the value of the lower envelope at any point is simply the minimum value of the elements in in the region defined by . It should be noted that the MAX and MIN filters produce new 2D matrices for upper and lower envelope surfaces from the given 2D data matrix, it does not alter the actual 2D data. Since smooth continuous surfaces for upper and lower envelopes are preferable, an averaging smoothing operation is carried out on both and employing the same window size used for corresponding order statistics filters. This averaging smoothing operation may be expressed as below:(11) where is the square region of size centered at any point of or , is the window width of the averaging smoothing filter and . The operations in (11) are arithmetic mean filtering that smoothes local variations in data. From the smoothed envelopes and , the mean or average envelope is calculated as in the original BEMD method given in Section 2.

As previously mentioned, FABEMD differs from BEMD in the way of formulating the upper and lower envelopes, and , and in restricting the number of iterations for each BIMF to one. In fact, one iteration per BIMF in FABEMD produces similar or better results than can be achieved by BEMD with more than one iteration. On the other hand, scattered data interpolation itself is an iterative process that fits a surface over the scattered data points in multiple steps. Though upper and lower envelope formation in FABEMD requires three steps: window size determination, getting the MAX (MIN) filter output, and averaging smoothing, all these operations can be done very fast using efficient programming routines; and the time required is much less than is required in the interpolation-based envelope estimation.

3.2.3. Illustration of Upper and Lower Envelopes Estimation

For illustration purposes of the envelope formation in FABEMD, let us consider the 2D data of Figure 1(a) and corresponding local maxima and minima maps of Figures 1(b) and 1(c). Window width obtained using Type-4 OSFW is 3. So, taking a window for MAX and MIN filters and applying them to the data matrix of Figure 1(a) results in the upper and lower envelope matrices given in Figures 3(a) and 3(b), respectively.

Figure 3: For data matrix of Figure 1(a): (a) upper envelope matrix using FABEMD before smoothing, (b) lower envelope matrix using FABEMD before smoothing.

The application of averaging smoothing operations to Figures 3(a) and 3(b) results in the smoothed upper and lower envelope matrices shown in Figures 4(a) and 4(b), respectively. The mean envelope matrix produced by averaging the matrices of Figures 4(a) and 4(b) is shown in Figure 4(c). For comparison purpose, corresponding matrices for UE, LE, and ME derived by using thin-plate spline surface interpolation to the maxima and minima maps are shown in Figure 5. Comparison of data in Figures 1(a), 4(c), and 5(c) reveals that the mean envelope derived by FABEMD method more closely matches the local mean of the given data. Since local mean subtraction is essential for the BEMD or FABEMD process to yield nearly zero local mean BIMFs, the FABEMD achieves this goal in as few as one iteration. It is shown in the literature [1, 3] that IMF or BIMF properties are retained when local mean is defined as the local mean of the upper and lower envelopes, not just the usual local mean as might be obtained by averaging the data using a spatial averaging filter smaller than the original size of the data matrix. Nevertheless, zero local envelope mean that also induces zero local mean yields well-characterized BIMFs.

Figure 4: For data matrix of Figure 1(a): (a) upper envelope matrix using FABEMD after smoothing, (b) lower envelope matrix using FABEMD after smoothing, and (c) mean envelope matrix obtained by averaging the data in (a) and (b).

Figure 5: For data matrix of Figure 1(a): (a) upper envelope matrix using BEMD with thin-plate spline interpolation, (b) lower envelope matrix using BEMD with thin-plate spline interpolation, and (c) mean envelope matrix obtained by averaging the data in (a) and (b).

To visualize the envelope formation for FABEMD more explicitly, let us consider a 1D signal for simplicity, given in Figure 6, where local maxima points are indicated by “” and local minima points are indicated by “x” that are obtained using a neighboring window. The sorted array of AMAXDA for this signal appears as , while the sorted array of AMINDA appears as . Using these distance arrays, OSFW for Type-1, Type-2, Type-3, and Type-4 appears to be 73, 79, 107, and 109, respectively. Taking Type-4 OSFW (i.e., ) as the width of the MAX and MIN filters and applying them to the 1D signal of Figure 6 results in the UE, LE and ME shown in Figure 7(a). The corresponding envelopes after applying smoothing averaging filter of the same size are displayed in Figure 7(b), and the same envelopes created by applying cubic spline interpolation to the maxima and minima maps are given in Figure 7(c). Figure 7(c) indicates the possibility of incorrect interpolation at the boundary and thus causing improper ME. The top waveforms in Figures 8(a) and 8(b) are the original 1D signal given in Figure 6, whereas the bottom waveform in Figure 8(a) is the result of ME subtraction in FABEMD method and the bottom waveform in Figure 8(b) is the result of ME subtraction in BEMD method. This illustration, along with the previous analyses, demonstrates the effectiveness of the proposed FABEMD method for BIMF or BEMC extraction.

Figure 6: A 1D signal and its local maxima and minima points.

Figure 7: (a) Envelopes using the proposed approach before smoothing, (b) envelopes using the proposed approach after smoothing, (c) envelopes using cubic spline interpolation.

Figure 8: (a) Original signal (top) and mean envelope subtracted signal (bottom) using FABEMD algorithm, (b) original signal (top) and mean envelope subtracted signal (bottom) using BEMD algorithm.

4. Simulation Results

The effectiveness of the FABEMD is investigated by implementing the algorithm for analyzing various images. The decomposed BEMCs resulting from FABEMD are compared with the BEMCs acquired using BEMD. Simulation results are reported for FABEMD with OSFW of Type-1 and Type-4. Although only one iteration for each BIMF is suggested in the FABEMD method, some results are also shown for more than one iteration to justify the adequacy of performing one iteration. Since the window sizes for order-statistics and smoothing filters are determined from the source image information, these sizes remain the same for all the iterations for the corresponding BIMF. FABEMD results are compared with the BEMD results obtained by thin-plate spline (TPS) interpolator, a radial basis function (RBF) that has been established as a good choice for BEMD [3–5]. For BEMD with RBF-TPS, SD criterion is employed as the fundamental stopping criteria with a threshold of 0.01, while the maximum number of allowable iterations (MNAI) is applied as additional stopping criterion to prevent over sifting [6, 15]. Additionally, in some cases BEMD results are also examined and reported for one iteration, to compare with the results of FABEMD with one iteration. To further limit the number of iterations and thus prevent over sifting in BEMD, SD defined by (1b) is considered for the simulation. It should be noted that the definition of SD affects the number of required iterations to achieve a given threshold and thus, the amount of sifting per BIMF; it does not have any contribution to the calculation of UE, LE, or ME in a particular iteration. Even though the complete space-spatial-frequency analysis using BHHT is investigated, the results are not shown in this paper. However, it is obvious that a good set of BIMFs will yield a good BHHT-based image representation. In the simulation, the maximum image size is limited to -pixel. Although FABEMD is capable of decomposing images of any size or resolution very fast (e.g., in few seconds or few minutes), BEMD is unable to do so. Since FABEMD results are compared with BEMD results for the same images, -pixel images help perform the task conveniently.

4.1. Analysis with Synthetic Texture Image

A synthetic texture image (STI) of -pixel size is taken, which is composed by adding three different components of the same size. For convenience of synthesizing, each synthetic texture component (STC) is generated from horizontal and vertical sinusoidal waveforms having different but closely spaced frequencies. The first STC consists of higher frequencies, the second STC consists of medium frequencies, and the last STC consists of very low frequencies. The STI and STCs are shown in Figure 9, while the diagonal intensity profiles of the STI and STCs are presented in Figure 10. Even if the addition of arbitrarily developed STCs in Figures 9(a) to 9(c) yields the original STI of Figure 9(d), application of BEMD or FABEMD to the STI of Figure 9(d) may not necessarily regenerate the STCs of Figures 9(a) to 9(c) (e.g., BEMC-1 may not be the same as STC-1), a consequence that can be attributed to the property of BEMD/FABEMD. Still, analysis of BEMD/FABEMD employing this synthetic texture provides a good performance indication of the algorithm. The OI among the original STCs and the global mean of each component are given in Table 1, which facilitates the comparison of the extracted BEMCs with the actual STCs. Since STC-1 and STC-2 are nearly symmetric with bipolar gray level values, their global mean should be close to zero as seen from Table 1. Thus, for STC-1 and STC-2 zero local envelope mean also implies zero global mean or vice versa.

Table 1: Global mean of STCs and their OI.

Figure 9: (a) Component 1 (STC-1), (b) component 2 (STC-2), (c) component 3 (STC-3), (d) original synthetic texture image (STI) obtained from addition of (a) to (c).

Figure 10: 1D diagonal intensity profiles of (a) STC-1, (b) STC-2, (c) STC-3, (d) STI.

Before demonstrating the final results of STI decomposition using FABEMD and BEMD, let us investigate the UE, LE, and ME of the STI generated by using different approaches. Figures 11(a) to 11(c) display the combined three-dimensional (3D) mesh plots of -pixel regions taken from the same locations of the original -pixel STI and the envelopes obtained by FABEMD with Type-4 OSFW, FABEMD with Type-1 OSFW, and BEMD with RBF-TPS interpolation, respectively. Figure 11 manifests the effectiveness of the proposed scheme of envelope estimation, which can very well replace the interpolation-based envelope estimation. Computation time of mean envelope estimation for the -pixel STI is also given in the parenthesis of the corresponding caption of Figure 11. In this case, it is noticed that the envelope estimation takes much shorter time with FABEMD than with BEMD. In general, OSFW increases for the later source images and hence the corresponding computation times of the envelopes also increase for the later BIMFs in the FABEMD process. On the other hand, the number of extrema points decreases for the later source images or ITS-BIMFs and therefore envelope estimation time decreases for later BIMFs in the BEMD process. The overall computation time in the BEMD process still remains much higher due to the iterative surface-fitting problem from the scattered data.

Figure 11: Mesh plots of -pixel regions taken from the -pixel STI and its UE, LE, and ME employing (a) FABEMD Type-4 OSFW (4.298072 seconds), (b) FABEMD Type-1 OSFW (3.717336 seconds), (c) BEMD RBF-TPS (193.124406 seconds).

Decomposition of the STI in Figure 9(d) is first conducted by applying FABEMD having Type-4 OSFW with . The resulting BEMCs and the summation of the BEMCs are displayed in Figure 12(a); and the corresponding diagonal intensity profiles are displayed in Figure 12(b). Figure 12 reveals the similarity of the BEMCs with the original STCs very well.

Figure 12: Decomposition of the STI using FABEMD with Type-4 OSFW () (a) BEMC-1 to BEMC-3 and summation of the BEMCs, (b) diagonal intensity profiles of BEMC-1 to BEMC-3 and summation of the BEMCs.

As mentioned in Section 2.2, in any approach of BEMD or FABEMD, the summation of the BEMCs will always return the original image back, except for the truncation/rounding error introduced at various steps of the process. This fact can be well verified from comparison of the STI and the summation of BEMCs in Figures 9, 10, and 12. Hence, showing the summation of BEMCs is excluded in the subsequent analyses of this paper.

The BEMCs of the STI obtained by applying FABEMD with Type-4 OSFW are displayed in Figure 13(a) for . The diagonal intensity profiles of the corresponding images in Figure 13(a) are shown in Figure 13(b). Because of the increased iterations, the stopping point SD for each BIMF decreases. This helps in attaining a first BIMF (BIMF-1), which is more similar to the original STC-1. But, due to the over sifting, an additional component appears that does not have any similarity to any of the original STCs. By looking at BEMC-3 of this decomposition, it may be inferred that this type of component may not have any significance in actual image processing applications. In fact, BEMC-3 and BEMC-4 may be combined to get a component similar to STC-3, although BEMC-3 contains some higher spatial scales compared to STC-3. Since the characteristics of diagonal intensity profiles for various BEMCs of the STI are now realized, displaying these profiles will be left out in the subsequent analyses.

Figure 13: Decomposition of the STI using FABEMD with Type-4 OSFW () (a) BEMCs (b) diagonal intensity profiles of BEMCs.

As a third example, the decomposed BEMCs employing FABEMD with Type-1 OSFW for are shown in Figure 14. Because Type-1 OSFW gives the minimum possible width from the distance matrix, it causes an increased level of sifting and thus a greater number of BIMFs/BEMCs (e.g., six BEMCs in this case). This reveals the fact that the selection of OSFW type can be made based on the image properties and desired applications.

Figure 14: BEMCs of the STI obtained by FABEMD with Type-1 OSFW ().

Application of FABEMD with Type-1 OSFW for generates seven BEMCs, which are displayed in Figure 15. This decomposition also shows the effect of over sifting and thus extraction of improper BEMCs in FABEMD. Due to the property of order statistics filter-based envelope estimation followed by a smoothing operation, over sifting in FABEMD may cause improper BEMCs as well.

Figure 15: BEMCs of the STI obtained by FABEMD with Type-1 OSFW ().

The STI is next decomposed using BEMD with RBF-TPS interpolation to compare the results with the new FABEMD method. The BEMCs resulting with are given in Figure 16 and the BEMCs resulting with are given in Figure 17. In both cases four BEMCs are generated, where BEMC-3 and BEMC-4 together may become similar to the STC-3. Due to the increased number of iterations, BEMC-1 appears better in Figure 17 than in Figure 16. In general, for BEMD, more than one iteration may generate more accurate BEMCs [3–5], although it is also better to limit the number of iterations without satisfying SD threshold criteria to prevent over sifting [6, 15, 21].

Figure 16: BEMCs of the STI obtained by BEMD with RBF-TPS ().

Figure 17: BEMCs of the STI obtained by BEMD with RBF-TPS ().

For additional performance assessment of the above six modes of the FABEMD/BEMD algorithm, the decomposition data of the STI are presented in Tables 2 to 5. Table 2 displays the number of obtained BEMCs, time taken for each algorithm, and OI for each case. In terms of time taken and orthogonality index, FABEMD employing Type-4 OSFW with appears to be the best choice for decomposing the STI considered in this paper. From Table 3 it is observed that one iteration of the process for each BIMF results in a comparatively higher stopping point SD. While higher SD implies nonzero local envelope mean for the corresponding BIMF, in the spatial scale-based decomposition, strict zero local envelope mean of a BIMF is not essential, which has been already established in the literature [6, 15]. This relaxation in turn allows prevention of over sifting by reducing the number of iterations, or not having a very low SD at the termination of the iterations. For FABEMD, increased iterations may cause erroneous envelopes and thus improper BIMFs, as mentioned previously. Table 4 represents the global mean of the BEMCs for all the FABEMD/BEMD modes applied to the STI. Like the first two symmetric and bipolar original STCs, the BIMFs are also expected to be symmetric and bipolar. For this reason, zero global mean of a BIMF should also indicate zero local mean or zero local envelope mean of it. In that sense, except the residue, other BEMCs (i.e., BIMFs) having nearly zero global mean can be considered good BIMFs. Hence, the method that produces a BIMF with higher global mean from the considered STI can be treated as poor. These features again designate FABEMD employing Type-4 OSFW with as a good choice for decomposition of the STI of Figure 9(d). Although the number of extrema in the source images for BEMD or FABEMD and the OSFW for each BIMF in FABEMD does not indicate any performance measure, these statistics are given in Tables 5 and 6 to provide some additional details of the processes.

Table 2: Comparison among various FABEMD/BEMD for the STI in terms of total number of BEMCs, total time required, and orthogonality index (OI).

Table 3: Comparison among various FABEMD/BEMD for the STI in terms of achieved stopping point SD for each BIMF.

Table 4: Comparison among various FABEMD/BEMD for the STI in terms of global mean of the BEMCs.

Table 5: Comparison among various FABEMD/BEMD for the STI in terms of number of extrema points in the source images for the corresponding BIMF, and the residue.

Table 6: Comparison among various FABEMD/BEMD for the STI in terms of order statistics filter width (OSFW).

4.2. Analysis with Real Images

Three real images are analyzed, in this section, to further investigate and compare the performance of FABEMD and BEMD. The first image is a -pixel region of a real texture image, D18, taken from the Brodatz texture set and shown in Figure 18(a) [24]. The second image is a subsampled -pixel Elaine image shown in Figure 18(b), while the third image is a -pixel noisy aurora image shown in Figure 18(c). The BEMCs generated from the D18 image, Elaine image, and aurora image by applying FABEMD with Type-1 OSFW and are shown in Figures 19, 21, and 23, respectively. On the other hand, the BEMCs generated from these same images by applying BEMD with RBF-TPS interpolation and 0 are shown in Figures 20, 22, and 24, respectively.

Figure 18: (a) -pixel region of Brodatz texture D18 [24], (b) -pixel Elaine image, (c) -pixel noisy aurora image.

Figure 19: BEMCs of D18 obtained by FABEMD with Type-1 OSFW ().

Figure 20: BEMCs of D18 image obtained by BEMD with RBF-TPS (0).

Figure 21: BEMCs of Elaine image obtained by FABEMD with Type-1 OSFW ().

Figure 22: BEMCs of Elaine image obtained by BEMD with RBF-TPS (0).

Figure 23: BEMCs of noisy aurora image obtained by FABEMD with Type-1 OSFW ().

Figure 24: BEMCs of noisy aurora image obtained by BEMD with RBF-TPS (0).

Because there are no ground-truth BEMCs of the considered real images, intuitive analysis using visual assessment is reported as the fundamental performance criterion in this paper. It is obvious from the BEMCs of real images that the FABEMD yields very well-defined BEMCs, which represent the image features at various spatial scales similar to, or better than, the BEMCs obtained from the BEMD method. Unwanted distortion and other artifacts may accompany the BEMCs when obtained via BEMD, which is apparent from the figures of BEMCs obtained by BEMD and given in Figures 20, 22, and 24; and thus they may not appear to be suitable for further image processing tasks. Although further evaluation of the texture decomposition can be reported by showing it for true texture analysis (e.g., texture classification, texture segmentation), achievement in having less distortion in the BEMCs of the texture image obtained with FABEMD is clearly visible in Figure 19. On the other hand, improvement of the BEMC quality for Elaine image and aurora image is obvious with FABEMD. For example, BEMCs of the Elaine image have no or less distortion, and more clearly reveal the edges and other characteristic features at different scales compared to the BEMCs obtained by BEMD. Similarly, observation of Figures 23 and 24 reveals that the noise is better separated into the first or first few BIMFs in FABEMD method than in BEMD method. This in turn should facilitate more efficient denoising of the aurora image using FABEMD compared to using BEMD. Since the obtained BEMCs are better in FABEMD, the complete BHHT analyses of images employing those BEMCs will be more effective. Preliminary studies on edge detection and noise removal using FABEMD show promising and significantly better performance compared to the analysis using BEMD, besides showing a dramatic improvement in the computation time. Since the objective of this paper is to provide the details of the FABEMD algorithm and its features, and to propose the algorithm for all types of applications, wherever BEMD-type processing may be used, the specific application-wise performance is not reported here.

Envelope estimation in FABEMD, employing order-statistics filters, is nearly independent of the image or texture pattern in terms of complexity and processing time; and the envelopes closely follow the image. But envelope estimation in the BEMD method, employing surface interpolation, is highly dependent on the maxima or minima maps, while the envelopes are not guaranteed to follow the image. In some cases when there are very few points in the maxima or minima maps, BEMD is prone to generate an erroneous surface and thus erroneous BEMCs. On the other hand, FABEMD is inherently free of boundary effects and overshoot-undershoot problems, and thus it does not require additional boundary processing. In Section 4.1, the time required for BEMD-based decomposition has been found to be higher than the time required for FABEMD-based decomposition of a simple and uniform STI. But for real images, the time taken by BEMD is even much higher than that required by FABEMD, which has been experienced in the example simulations presented in this section as well. While FABEMD takes only a few minutes, BEMD takes many hours, even for a very few iterations performed per BIMF. This problem hinders the application of BEMD in many practical cases. Adaptability achievable through the selection of OSFW is another supplementary feature of the FABEMD process. This feature allows various possible sets of BEMCs from the same image, which in turn helps in optimizing the image processing need by providing the opportunity of selecting an appropriate set. Even though this type of adaptability is also available from BEMD by means of selecting different types of interpolation, the associated problems of interpolation may not render the utilization of this method successfully.

Although various possibility of OSFW provides adaptability, sometimes it may impose some difficulty in applying the FABEMD algorithm; because, to achieve the desired decomposition, a trial and error selection procedure may have to be performed to find the required value for OSFW. On the other hand, the requirement for manipulation of for a BIMF, when the calculated value does not appear larger than the previous BIMF mode, may impose additional complexity. Hence, reducing the above-mentioned difficulties can be worked out in future. Also, applying the FABEMD algorithm for various real image processing tasks and reporting the methodologies and the corresponding results individually with comparison to the results employing other BEMD approaches will be interesting.

5. Conclusion

BEMD is a potential image processing algorithm. To boost increased application of this algorithm for image processing applications, a fast, time efficient, and effective method is essential. This fact motivated the formulation of a fast and adaptive BEMD, abbreviated as FABEMD, described in this paper. In FABEMD, the envelope estimation method of regular BEMD is modified by replacing the 2D surface interpolation by an order-statistics-based filtering followed by a smoothing operation. A number of window sizes can be selected for the order statistics and smoothing filters, all of which are data driven and thus making the process adaptive. The simple change in the envelope estimation procedure provides a tremendous enhancement of the algorithm in terms of computation time. The proposed FABEMD has been tested for decomposing various images, some of which have been reported in this paper. Simulation results demonstrate the usefulness of this novel FABEMD approach for BEMD-based image decomposition. FABEMD enables the decomposition of images with any dimensions in a very short period of time, while the application of BEMD is still limited to smaller images. Beside reducing the computation time, this novel approach also ensures a more accurate estimation of the BIMFs in some cases. It is believed that FABEMD can be a perfect alternative to the regular BEMD and will play a very significant role in this area.

References

1. N. E. Huang, Z. Shen, S. R. Long, et al., “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proceedings of the Royal Society A, vol. 454, no. 1971, pp. 903–995, 1998.
2. N. E. Huang, Z. Shen, and S. R. Long, “A new view of nonlinear water waves: the Hilbert spectrum,” Annual Review of Fluid Mechanics, vol. 31, pp. 417–457, 1999.
3. J. C. Nunes, Y. Bouaoune, E. Deléchelle, O. Niang, and Ph. Bunel, “Image analysis by bidimensional empirical mode decomposition,” Image and Vision Computing, vol. 21, no. 12, pp. 1019–1026, 2003.
4. J. C. Nunes, S. Guyot, and E. Deléchelle, “Texture analysis based on local analysis of the bidimensional empirical mode decomposition,” Machine Vision and Applications, vol. 16, no. 3, pp. 177–188, 2005.
5. A. Linderhed, “2-D empirical mode decompositions in the spirit of image compression,” in Wavelet and Independent Component Analysis Applications IX, vol. 4738 of Proceedings of SPIE, pp. 1–8, Orlando, Fla, USA, April 2002.
6. C. Damerval, S. Meignen, and V. Perrier, “A fast algorithm for bidimensional EMD,” IEEE Signal Processing Letters, vol. 12, no. 10, pp. 701–704, 2005.
7. Y. Xu, B. Liu, J. Liu, and S. Riemenschneider, “Two-dimensional empirical mode decomposition by finite elements,” Proceedings of the Royal Society A, vol. 462, no. 2074, pp. 3081–3096, 2006.
8. Z. Liu, H. Wang, and S. Peng, “Texture classification through directional empirical mode decomposition,” in Proceedings of the 17th IEEE International Conference on Pattern Recognition (ICPR '04), vol. 4, pp. 803–806, Cambridge, UK, August 2004.
9. Z. Liu, H. Wang, and S. Peng, “Texture segmentation using directional empirical mode decomposition,” in Proceedings of IEEE International Conference on Image Processing (ICIP '04), vol. 1, pp. 279–282, Singapore, October 2004.
10. S. R. Long, “Applications of HHT in image analysis,” in Hilbert-Huang Transform and Its Applications, N. E. Huang and S. S. P. Shen, Eds., World Scientific, River Edge, NJ, USA, 2005.
11. Z. Liu and S. Peng, “Boundary processing of bidimensional EMD using texture synthesis,” IEEE Signal Processing Letters, vol. 12, no. 1, pp. 33–36, 2005.
12. N. E. Huang, M.-L. C. Wu, S. R. Long, et al., “A confidence limit for the empirical mode decomposition and Hilbert spectral analysis,” Proceedings of the Royal Society A, vol. 459, no. 2037, pp. 2317–2345, 2003.
13. G. Rilling, P. Flandrin, and P. Gonçalves, “On empirical mode decomposition and its algorithms,” in Proceedings of IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing (NSIP '03), Grado, Italy, June 2003.
14. C. Junsheng, Y. Dejie, and Y. Yu, “Research on the intrinsic mode function (IMF) criterion in EMD method,” Mechanical Systems and Signal Processing, vol. 20, no. 4, pp. 817–824, 2006.
15. M. Shen, H. Tang, and B. Li, “The modified bidimensional empirical mode decomposition for image denoising,” in Proceedings of the 8th IEEE International Conference on Signal Processing (ICSP '06), vol. 4, Beijing, China, November 2007.
16. Y. Washizawa, T. Tanaka, D. P. Mandic, and A. Cichocki, “A flexible method for envelope estimation in empirical mode decomposition,” in Proceedings of the 10th International Conference on Knowledge-Based Intelligent Information and Engineering Systems (KES '06), B. Gabrys, R. J. Howlett, and L. C. Jain, Eds., pp. 1248–1255, Bournemouth, UK, October 2006.
17. R. Srinivasan, R. Rengaswamy, and R. Miller, “A modified empirical mode decomposition (EMD) process for oscillation characterization in control loops,” Control Engineering Practice, vol. 15, no. 9, pp. 1135–1148, 2007.
18. Z. K. Peng, P. W. Tse, and F. L. Chu, “An improved Hilbert-Huang transform and its application in vibration signal analysis,” Journal of Sound and Vibration, vol. 286, no. 1-2, pp. 187–205, 2005.
19. B. Xuan, Q. Xie, and S. Peng, “EMD sifting based on bandwidth,” IEEE Signal Processing Letters, vol. 14, no. 8, pp. 537–540, 2007.
20. Z. Liu, “A novel boundary extension approach for empirical mode decomposition,” in Proceedings of the International Conference on Intelligent Computing (ICIC '06), D. Huang, K. Li, and G. W. Irwin, Eds., pp. 299–304, Kunming, China, August 2006.
21. B. Shen, “Estimating the instantaneous frequencies of a multicomponent AM-FM image by bidimensional empirical mode decomposition,” in Proceedings of IEEE International Workshop on Intelligent Signal Processing (WISP '05), pp. 283–287, Faro, Portugal, September 2005.
22. S. Kizhner, K. Blank, T. Flatley, N. E. Huang, D. Petrick, and P. Hestnes, “On certain theoretical developments underlying the Hilbert-Huang transform,” in Proceedings of IEEE Aerospace Conference, p. 14, Big Sky, Mont, USA, March 2006.
23. R. C. Gonzalez and R. E. Woods, Digital Image Processing, Pearson Education, Upper Saddle River, NJ, USA, 2nd edition, 2002.
24. P. Brodatz, Textures: A Photographic Album for Artists and Designers, Dover, New York, NY, USA, 1966.