## Abstract

Empirical mode decomposition (EMD) is a powerful tool for analysis of non-stationary and nonlinear signals, and has drawn significant attention in various engineering application areas. This paper presents a finite element-based EMD method for two-dimensional data analysis. Specifically, we represent the local mean surface of the data, a key step in EMD, as a linear combination of a set of two-dimensional linear basis functions smoothed with bi-cubic spline interpolation. The coefficients of the basis functions in the linear combination are obtained from the local extrema of the data using a generalized low-pass filter. By taking advantage of the principle of finite-element analysis, we develop a fast algorithm for implementation of the EMD. The proposed method provides an effective approach to overcome several challenging difficulties in extending the original one-dimensional EMD to the two-dimensional EMD. Numerical experiments using both simulated and practical texture images show that the proposed method works well.

## 1. Introduction

Empirical mode decomposition (EMD) was originally developed for analysis of one-dimensional non-stationary and nonlinear signals (Huang *et al*. 1998). It decomposes a signal into a finite sum of intrinsic mode functions (IMF) that generally allow well-behaved Hilbert transforms. This makes it possible to construct a time–frequency representation, called the Hilbert spectrum, using instantaneous frequency. EMD combined with Hilbert spectrum, called the Hilbert–Huang transform (HHT), has many advantages over the traditional time–frequency analysis techniques in the adaptability, the capability of analysing nonlinear signals and the high time–frequency resolution. In recent years, HHT has been applied with great success in various application areas.

Beyond time–frequency analysis, EMD has a very important feature in that it is a fully adaptive multiscale decomposition (Wu & Huang 2004; Flandrin *et al*. 2004). This is because EMD operates on the local extremum sequence and the decomposition is carried out by direct extraction of the local energy associated with the intrinsic time-scales of the signal itself. This is different from the wavelet-based multiscale analysis that characterizes the scale of a signal event using pre-specified basis functions. Owing to this feature, EMD is highly promising in dealing with other problems of a multiscale nature. The present research is an attempt in this respect and aims at providing an EMD method for two-dimensional data analysis. Although this problem has received attention recently (Linderhed 2002; Nunes *et al*. 2003), there are several challenging difficulties that need to be overcome. One of them is the computation efficiency. For a medium-sized two-dimensional data, e.g. a 512×512 image, the number of the local extrema can be tens of thousands. Owing to the iterative nature of the EMD method, the decomposition of such a dataset is rather time-consuming and could be unacceptable for many applications. While converting a two-dimensional data to one dimension and then employing the one-dimensional EMD could be an efficient approach to deal with some image processing problems (Yang *et al*. 2004), development of a fast two-dimensional EMD is useful for more general applications.

In development of the two-dimensional EMD method, while the sifting process of the original EMD can be directly used, the key step is to generate the local mean surface of the two-dimensional data. The present research extends the basic idea of the recently developed B-spline EMD (Chen *et al*. 2006) and uses finite-element basis functions, which replace B-splines in the one-dimensional case, to construct the local mean surface of the data. Specifically, for a dataset defined on a two-dimensional domain, we first generate a triangular mesh on the domain with their vertices being the local extrema and the saddle points of the given data. Then, we construct a linear combination of the set of linear shape functions defined on the triangular mesh. The coefficients of the shape functions are obtained by applying a generalized low-pass filter to the local extremum and saddle point set. The local mean surface of the data is finally obtained by smoothing the linear combination with bi-cubic spline interpolation. This procedure can be implemented efficiently. In addition, the determination of the parameters requires less user involvement than the existing methods. It was brought to our attention by a referee that a fast algorithm for two-dimensional EMD was independently developed in a recent paper (Damerval *et al*. 2005), which was published while our paper was under review. The paper extends the one-dimensional EMD algorithm of Huang to the two-dimensional case by using cubic interpolation on triangles to generate the upper and lower envelope surfaces. The present research differs from the work in Damerval *et al*. (2005) mainly in that we generate the local mean surface of a two-dimensional dataset directly from the characteristic data points rather than from the upper and lower envelopes. This overcomes the problem of possible overshootings between the upper and lower envelopes. Our method avoids constructing two different two-dimensional interpolating surfaces (the upper and lower envelopes), which is normally a difficult task and requires much computational cost. In addition, the characteristic data points in our method include not only the local maxima and minima, but also the saddle points, which are a distinct feature of two-dimensional data. We apply the proposed method to some simulated and practical texture image data. The results indicate that it works well for our collection of test images.

The paper is organized into six sections. In §2, we present a general scheme of EMD that forms the data local mean as the moving average of non-oscillated basis functions. Based on this scheme, the two-dimensional EMD method is developed in §3. In §4, we describe an algorithm for fast evaluating values of the local mean surface in the two-dimensional EMD. Then, in §5, we apply the proposed method to analyse simulated two-dimensional datasets and practical images. We summarize the research in §6.

## 2. A general setting for EMD

A crucial step in the EMD method for a signal is to represent its *local mean*. In the original EMD (Huang *et al*. 1998), the local mean of the data is defined as the mean of the upper and lower envelopes, which are computed by using the cubic spline interpolation. Using the mean of the upper and lower envelopes as a local mean is motivated by the requirement to compute instantaneous frequency using Hilbert transform. This is altered in Chen *et al*. (2006) for the convenience of studying the mathematical properties of EMD, where the local mean is represented as a moving average of B-splines with the knots being the local extrema of the data. It has been shown that the B-spline EMD has a comparable performance with the original EMD for all our test cases (Chen *et al*. 2006; Riemenschneider *et al*. 2005; Liu *et al*. 2006). In fact, the basic idea of Chen *et al*. (2006) is also applicable to high-dimensional cases. In this section, we present a general setting for EMD. This generalization serves two purposes. First, it helps us understand the intrinsic properties of the EMD method and second, it will allow us to develop a finite element-based two-dimensional EMD method.

We now describe the general setting of EMD. Let , for . We choose a compactly supported basis , , for a subspace of , which are adapted to the given data. These basis functions will be used to construct the local mean surface of the given data. Suppose that the setis the collection of characteristic points of *f*, which include the extreme points for the one-dimensional EMD, the extreme and saddle points for the two-dimensional EMD, as defined in §3. We define a smoothed functional of by(2.1)where *S* is a generalized low-pass filter and represents a generalized filtering operation. A class of generalized low-pass filters and high-pass filters having vanishing moments are found in Chen *et al*. (2002) for the purpose of developing collocation methods for solving integral equations. Basically, we specify the generalized low-pass filter as a set of positive-valued weights and the generalized filtering as a weighted sum of and its neighbouring local extrema. As such, the produced functional will be smoother than since the local variation of is averaged.

With the basis functions and the smoothed functionals of the local extrema, we define the local mean of the data *f* as(2.2)where are the basis functions. In fact, *m* is a quasi-interpolation to the original data *f*. Note that we avoid using the ‘upper envelope’ and ‘lower envelope’ in the definition of the local mean, because mathematically they are not well-defined by using spline interpolation.

The essence of the EMD method is to subtract the local mean from the data so as to decompose the data into a high-frequency and a low-frequency component (i.e. the local mean). The EMD method decomposes a signal into a finite sum of IMF. In the one-dimensional case, an IMF is originally defined in terms of the difference between the number of zero crossings and the number of extrema of the signal, so that the instantaneous frequency of an IMF is well defined by the Hilbert transform. In a general case, in particular in the two-dimensional case, we will not impose a specific definition of an IMF. It should be determined by a specific stopping condition in the sifting process for a specific application.

We now describe a general EMD algorithm. We extract the first IMF by the following steps.

Find the local extrema .

Compute the smoothed set of the local extrema using equation (2.1).

Compute the local mean

*m*using equation (2.2).Compute . If

*h*satisfies a given stopping condition, stop. Otherwise, treat*h*as the data and iterate on*h*.

We denote by the first IMF and specify as the first residue. The next IMF, , is then selected by applying the above procedure to the first residue . This process is repeated until a satisfactory result is obtained.

The procedure generates a sequence of IMFs, , and a residue function . The basic idea of the EMD is to decompose a signal into the sum of these IMFs and the residue function, so that catches the highest frequency of *f*, the second highest frequency of *f* and catches the lowest frequency of *f*. Specifically, this procedure yields a decomposition(2.3)

This general setting covers the envelope approach, the original EMD presented in Huang *et al*. (1998), the B-spline approach developed in Chen *et al*. (2006) and the two-dimensional EMD by using finite elements, which we will present in this paper. It allows us to think of the EMD in more general terms. A different way of defining the local mean gives a different method for construction of EMD. The local mean can be defined by interpolation, by quasi-interpolation, or by other approximation approaches.

For the original EMD developed in Huang *et al*. (1998), the basis functions are cubic splines determined by the extreme points of the signal. The smoothed functional is the coefficient of the cubic splines defined by the cubic spline interpolation. In this case, both basis functions and the smoothed set are determined implicitly.

In the B-spline EMD, the basis functions and the smoothed functional are defined explicitly. In fact, the basis functions are chosen as B-splines with the knots being the extreme points of the signal. The generalized low-pass filter is defined as the binomial sequence and the smoothed sequence of the local extrema is the moving average(2.4)where *k* is the order of the B-spline function , and , are the extreme points inside the support of the B-spline . The local mean defined in this way has certain advantages over the approaches using envelopes. It overcomes the overshooting problem of the upper and lower envelopes (see Chen *et al*. 2006).

In the two-dimensional EMD that we will describe in §3, the basis functions are smoothed linear finite-element shape functions and the generalized low-pass filter is chosen as a weighted average of the function values at the extreme points around the point of interest. This construction is a natural extension of the one-dimensional B-spline EMD to the two-dimensional case. It also overcomes the same overshooting problem if the upper and lower envelopes of interpolation are used to construct the local mean surface.

For the one-dimensional B-spline EMD and the two-dimensional finite-element EMD, the computation cost is less than the envelope-based approach, since it requires evaluating only once the function value at a point in each iteration of the sifting process. Note that while the determination of the local mean for the two-dimensional data is different from the one-dimensional case, we use the same sifting procedure to decompose the data and obtain a decomposition (2.3).

In Liu *et al*. (2006), we have indicated that for the one-dimensional B-spline case, the local mean, defined by equations (2.2) and (2.4), represents the low-frequency part of a signal. In general, since the setis a smoothed version ofthe local mean will be a smoothed version of the data *f* if the basis functions themselves do not introduce more oscillations. In the two-dimensional case presented in §3, the basis functions are selected as a *smoothed* linear shape functions to ensure that no new oscillations are introduced. Hence, by smoothing the edges of the shape functions, the produced local mean surface also represents the low-frequency part of the data. The subtraction of the low-frequency part *m* from the data *f* yields the high-frequency part. In the next iteration of the sifting process, a new low-frequency part is generated from the current high-frequency part. When the stopping condition is satisfied, we obtain a high-frequency part, i.e. the first IMF, and the residue, which is actually the sum of the low-frequency parts generated in all the iterations of the sifting process. The next IMF is extracted in the same way but from the residue. When all the IMFs are obtained, the data will be decomposed into a number of frequency bands with the low-frequency bands corresponding to the large scales, while the high-frequency bands correspond to the small scales.

## 3. The finite-element method for the two-dimensional EMD

In this section, we follow the general setting for EMD described in §2 to propose a two-dimensional EMD based on the finite elements. A key step is to construct a local mean surface using the *characteristic points* (which we will define later) of the original function (or data). There are two major building blocks in this construction, which consists of the finite-element basis functions and the smoothed data. Our idea is to use the characteristic points of the original data to define a triangle partition of the domain and to build finite-element basis functions on the partition. The finite-element basis functions are the analogy of the B-spline functions used in the one-dimensional EMD developed in Chen *et al*. (2006). We propose to use a local weighted average of the values of the function at the characteristic points near the point under consideration as the smoothed data. This construction is a natural two-dimensional extension of the one-dimensional B-spline based EMD.

The first crucial issue that we face is what kind of features of the original function is needed to construct the local mean. It is well known that in one-dimensional EMD, the local maxima and local minima are important. In the two-dimensional case, in addition to the local maxima and local minima in the usual sense, we also need to consider the saddle points, since they are simultaneously a local maximum and local minimum point evaluated in different orientations, and they also give important features about the local variation of the original function. For this reason, we do not distinguish the three types of points and instead we treat them in the same way. A point is called a *characteristic point* if it is either a local maxima, local minima or a saddle point. As such, the decomposition will assign reasonably the local oscillated components characterized by saddle points to the high-frequency part, i.e. the IMF part, rather than the local mean part.

We now describe our idea in detail. Let be a closed polygonal domain and . Suppose thatis the collection of the characteristic points of the function *f* on , where is an appropriate index set. As illustrated in figure 1, we partition the domain into a triangular mesh with vertices being the characteristic points. This can be realized, for example, by using the Delaunay method (Sapidis & Perucchio 1991). In implementation, we also include the corner points of the domain in the set of characteristic points. This ensures that the triangular mesh covers the whole domain . In such a triangular mesh, any triangle does not overlap with any other triangles in the mesh. Moreover, a vertex of a triangle is not in the interior of an edge of another triangle in the mesh. Hence, if is not on the boundary of , then there is a finite number of points , such that they are the vertices of the polygon surrounding and no other points of , except located interior to the polygon. The set is called the set of neighbouring characteristic points around . For , we let be a triangle with the vertex and two other points in and let be the cardinality of . We denote the polygon around by and observe thatWhen is on the boundary of , one can extend the domain , so that is an interior point of the extended domain. One approach for the extension is the symmetric extension. As shown in figure 1, after the partition, any characteristic point on the boundary of is a vertex of the triangles around . The symmetric extension is done by mirroring with respect to the boundary all the triangles that have one vertex at . Specifically, due to the simplicity of the linear finite element, in our actual implementation, we do not need to perform the extension because the polygon around obtained from the symmetric extension is symmetric with respect to the boundary and the subsequent computation can be done only on the side within the domain . We use this simple approach in the numerical examples presented in §5. However, we should point out that handling the boundary condition is still a challenging issue for EMD and has also received attention in the study of its two-dimensional extensions (e.g. Damerval *et al*. 2005). While currently a perfect extension method is generally unavailable, the corresponding computation cost is a factor preferably considered in designing an extension approach.

Associated with the triangle partition, we build our finite-element basis functions. Specifically, we assign to each , a basis shape function . Outside the polygon , is equal to zero and on each it is a linear polynomial satisfying the interpolation conditions(3.1)Note that the function and it is a piecewise linear polynomial defined on having the support being the polygon . Such functions are finite-element shape functions and they have been used extensively in numerical solutions of partial differential equations. One can also construct smoother (such as or ) shape functions using higher-order polynomials. For more about such constructions, properties of finite-element functions and their use in numerical solutions of boundary-value problems, the reader is referred to Ciarlet (1976).

We next describe a construction of the smoothed dataset. We choose the smoothed data as a weighted average of the values of the signal *f* at the neighbouring characteristic points of the point . For example, we may specifically choose(3.2)The parameter *α* controls the degree of the smoothing and it is chosen empirically. In the numerical examples presented in this paper, it is specified as . The formulation (3.2) is called the *Laplacian operator* in the literature (cf. Taubin 1995; Kobbelt *et al*. 1998). One can also use other forms of local weighted average of the values of the signal *f* at the characteristic points around the point to define the smoothed data .

According to equation (2.2), we generate a mean surface(3.3)Let *T* be a triangle element of the triangularization of domain and suppose that its vertices are . Only three shape functions , , are non-zero at an interior point of *T*. Hence, for an interior point *p* of *T*, we have that(3.4)Clearly from this formula, the mean can be evaluated locally. In the two-dimensional EMD, a crucial issue is computation efficiency, since it is necessary to process a large number of extreme points even for a medium-sized two-dimensional dataset. Owing to the construction of (3.4), the evaluation of the function values of , which represents the major computation cost of the two-dimensional EMD, can be efficiently implemented.

A *smooth* local mean *m* is desirable for application purposes to avoid introducing additional new oscillations. The surface generated by equation (3.3) or (3.4) is continuous but not smooth (i.e. it is in but not in ), because linear shape functions are used as basis functions. This may consequently introduce additional sharp structures into the data. When a smooth mean surface is needed, we may choose to use smooth higher-order finite elements to overcome this problem. Since using higher-order finite elements may result in a significant increase in the computational cost, an alternative idea is to apply a smoothing filter to in order to obtain a new smooth mean *m*, avoiding the use of higher-order finite elements. Specifically, we will use a specially designed bi-cubic spline interpolation *m* of , so that (i) the new mean *m* is smooth, (ii) it preserves the crucial properties of and (iii) the additional computational cost is as small as possible.

We now describe a construction of the interpolation. Consider a two-dimensional dataset sampled with a rectangular sampling grid. First, we compute the number of the characteristic points on each row, denoted by , , where is the total number of rows of the grid. Let and be the maximum and minimum value of , , respectively. Then we pick up a row if its number of the characteristic points is greater than given by(3.5)where is a parameter specified empirically between 3 and 6. The first and last row of are always selected no matter whether their numbers of characteristic points are greater or less than . The columns of are processed in the same way.

The above procedure will produce a non-uniform and relatively sparse grid on the domain of the data. The knots for the spline interpolation are then determined as the characteristic points on each of the selected rows and columns, as well as the intersection points of the selected rows and columns. For the knots which are not characteristic points, we use the fast algorithm described in §4 to evaluate the values of . Our local mean surface *m* is finally specified as the smoothed function obtained from the bi-cubic interpolation.

Now we summarize the method of the two-dimensional EMD. Suppose that we have a dataset *f*. The first IMF of *f* is obtained as follows.

Find the local extrema and saddle points of

*f*.Form a triangular mesh using the Delaunay method.

Smooth the characteristic point set using the Laplacian operator of equation (3.2).

Pick up the rows and columns on which the number of the characteristic points is greater than the value specified in equation (3.5), determine the knots for the bi-cubic interpolation and generate the local mean surface .

Compute .

If(3.6)where is the sifting result in the

*k*th iteration, and*SD*is typically set between 0.2 and 0.3, stop and we obtain an IMF. Otherwise, treat*h*as the data and iterate on*h*through steps (i)–(vi).

Denote by the first IMF and setwhich is the first residue. The algorithm proceeds to select the next IMF by applying the above procedure to the first residue . This process is repeated until the last residue has at most one extremum or becomes constant. The original data can then be represented as(3.7)Equation (3.7) gives a two-dimensional EMD based on the finite elements.

The computation complexity in step (ii) of the algorithm is (Watson 1981), where *m* is the number of characteristic points of the current residue in a sifting iteration. Depending on the size and composition of the analysed image, this number may be large for the first few IMFs, but reduces rapidly for the remaining IMFs.

## 4. Fast evaluation of the local mean surface values

The method described in §3 heavily depends on an efficient evaluation of the local mean surface. We propose in this section a fast algorithm with linear computational complexity for evaluation of the local mean surface. Specifically, when evaluating the values of in equation (3.3), the most time-consuming procedure is to determine which triangle in the mesh contains the point *p* if *p* is not a node. Here, we present a fast algorithm to do this based on the finite-element analysis technique described in Kobbelt *et al*. (1998). The complexity of this algorithm is linear, i.e. , where *N* is the total number of the data points.

As shown in figure 2, we denote by *D* the domain of the original data *f*. It is divided into equally spaced rows and columns with a unit distance between the neighbouring columns and rows. For each triangle element *T* in the mesh, we first find the smallest rectangle that contains *T* according to the positions of the three nodes of *T*. Then, for each point *p* in , we determine whether it is in *T* as follows.

As shown in figure 2, the point *p* and every two nodes of *T* form a sub-triangle. Denote these sub-triangles by , . Let the coordinates of the three nodes of *T* be , and the coordinate of the point *p* be . Then, the direct areas of , and *T* can be computed byIt can be easily shown that the sufficient and necessary condition for a point *p* to be contained in the triangle *T* is(4.1)For , we introduce the notations(4.2)Then, we have a condition equivalent to (4.1) that has the form(4.3)It indicates that if the condition (4.3) is satisfied, then the point *p* is included in *T*, and vise versa.

Since the union of all the triangles covers *D*, the union of all the rectangles then includes *D*. Therefore, for any point in *D*, the triangle it belongs to can be determined using the above process. One can easily show that the complexity of the computation is linear.

## 5. Application examples

In this section, we present several examples of texture analysis to demonstrate the performance of the proposed two-dimensional EMD.

We first analyse a texture image, shown in figure 3*a*, taken from the Brodatz texture album (Brodatz 1966). This image appears to contain pattern structures of different scales. Our purpose is to separate these pattern structures, which is of interest in texture understanding and classification applications. In the analysis, the image is decomposed into three IMFs and a residue, as shown in figure 3*b–d*. For the two-dimensional EMD, it is not always needed to decompose the data to the final level, where the residue contains at most one local extremum, if the Hilbert transform is not performed. In the present example, the residue is not further decomposed, as it characterizes better the pattern structure of the largest scale. The three IMFs in figure 3*b–d* and the residue in figure 3*e* represent the pattern structures from finer to courser scales.

To further demonstrate the capability of the proposed method in multiscale texture decomposition, we analyse a simulated image for which we know the composition *a priori*. The image is shown in figure 4*a*. It consists of three components of different scales shown in figure 4*b–d*. Figure 4*e–g* shows the two IMFs and the residue obtained using the proposed method, corresponding to the small, medium and large-scale components, respectively. Comparing the third row in figure 4 with the second row, one can see that the proposed method has properly separated the components of different scales.

In the above simulated example, if we reconstruct an image using only the IMFs, we can then obtain an image without containing the tendency component, as shown in figure 5. This suggests an application of the two-dimensional EMD for removing inhomogeneous illumination, as also discussed in Nunes *et al*. (2003). In fact, the component produced by inhomogeneous illumination can often be modelled as a low-order polynomial without oscillation (Rowley 1999). It is thus contained in the residue of the EMD decomposition, and therefore can be removed through the reconstruction. Such an application is of interest in image recognition.

Now we present a more practical application example of the proposed method. It is the detection of the raw textile defects occurring in textile manufacturing. This is again a multiscale decomposition problem. Since in many cases the defects and the textile textures have different scales, we can thus use the two-dimensional EMD to separate these textile textures. This will then enhance the defects in the image and make them more noticeable. Figure 6*a* shows a textile image containing a defect called thick bar (Murino *et al*. 2004). The defect is rather obscure and could be neglected in the detection. The image is decomposed into four IMFs and one residue. In figure 6*b*, we present the reconstructed image of the IMFs, which represents mainly the textile texture. The residue image is presented in figure 6*c*. It shows clearly the existence of the defect. Figure 6*d* shows another image, where the defect is a thin bar (Murino *et al*. 2004). The image is processed in the same way as the one in figure 6*a*. Again, the defect has been enhanced successfully in the residue image, as shown in figure 6*f*.

## 6. Conclusion

In this paper, we have developed an EMD method for two-dimensional data analysis. This method is based on our general EMD scheme that represents the local mean of the data as a linear combination of non-oscillated basis functions. Specifically, the basis functions here are the two-dimensional linear functions smoothed with bi-cubic spline interpolation. The coefficients of the basis functions in the linear combination are the local extrema of the data smoothed using a generalized low-pass filter. Our method develops a fast algorithm based on the principle of finite-element analysis. We have applied the proposed method to simulated and practical texture image multiscale decomposition and raw textile defect detection. The results show that the two-dimensional EMD is particularly suitable to deal with texture images. It is anticipated that besides texture images the proposed method could also be useful to process other two-dimensional data with similar structures. However, as the principle of the EMD is based mainly on the local extrema of the data, it is thus not sensitive to edge features. For the images where edges are the major features, further research is still needed.

## Acknowledgments

The authors would like to thank Prof. Vittorio Murino of the Department of Computer Science, University of Verona, Italy, for providing the textile images used in this paper. Supported in part by National Science Foundation under grants 0314742 and 0312113, by National Aeronautics and Space Administration under grants NAG5-5364 and NGT4-52431, and by the Department of Energy under grant DE-FG02-04ER46116. Y.X. was supported in part by the US National Science Foundation under grant CCR-0407476, by the Natural Science Foundation of China under grant 10371122, and by the Education Ministry of the People's Republic of China under the Changjiang Scholar Chair Professorship Program through Zhongshan (Sun Yat-sen) University.

## Footnotes

- Received September 20, 2005.
- Accepted February 21, 2006.

- © 2006 The Royal Society