Abstract: This study presents a novel bi-dimensional EMD (Empirical mode decomposition) algorithm. It firstly transforms the given 2D data (e.g., an image) into a function defined on a triangular mesh and then uses a novel interpolation method based on bi-Laplace operator to generate the upper and lower envelopes in the sifting process. Essentially, this method still uses thin plate spline interpolation to generate envelopes but it can be implemented much faster than the traditional method based on thin plate spline interpolation because only a sparse linear system needs to be solved when generating an envelope in the sifting process. Furthermore, the proposed bi-dimensional EMD algorithm can be applied in image enhancement and the application examples prove that the image enhancement algorithm based on this bi-dimensional EMD approach is efficient and better than some classical image enhancement algorithms in contrast and lightness controlling.
INTRODUCTION
EMD (Empirical mode decomposition) is developed by Huang et al. (1998) and Huang and Wu (2008) which can adaptively decompose a 1-dimensional (1D) nonlinear and non-stationary signal into a finite sum of different frequency components. It is very important in signal processing and has been found useful applications in many areas such as signal denoising, fault detection, geophysical studies, atmospheric and climate studies, oceanographic studies and so on (Huang et al., 1998; Huang and Wu, 2008).
This study focuses on decomposing 2-dimensional (2D) data (e.g., images) by EMD technique, i.e., bi-dimensional EMD (BEMD). Currently, several methods about bi-dimensional EMD have been proposed. Nunes et al. (2003, 2005) generated envelope surfaces by radial basis function interpolation and used the Riesz transform for image analysis. Linderhed (2005) used a thin plate spline (a special case of radial basis function) for surface interpolation to develop bi-dimensional EMD algorithm for an image compression scheme. However, the determination of the thin plate spline is computationally expensive because a full linear system with the same number unknowns as data points must be solved. Damerval et al. (2005) proposed a way based on Delaunay triangulation and piecewise cubic polynomial interpolation to obtain envelope surfaces. Xu et al. (2006) used a mesh fitting method based on finite elements to represent the local mean surface of the data. Although, the methods based on Delaunay triangulation and finite-element interpolation can decompose the original signal much faster than the others, they are still time consuming (Bhuiyan et al., 2009).
In addition, each row and/or each column of 2D data can be processed, respectively by 1D EMD which makes it a faster process (Han et al., 2002). Unfortunately, this 1D implementation may generate uncorrected bi-dimensional intrinsic mode function components as it ignores the correlation among the rows and columns of 2D data.
Accounting for thin plate spline is a natural extension of cubic spline used in 1D EMD interpolation to 2D, this study presents a novel bi-dimensional EMD algorithm using thin plate spline interpolation based on bi-Laplace operator. Compared to the traditional approach based on thin plate interpolation (Linderhed, 2005), this approach can be implemented much faster because only a sparse linear system needs to be solved when generating an envelope in the sifting process. Experiment results and the application examples in image enhancement show that the proposed method is robust and efficient.
1D EMD
The 1D EMD algorithm proposed by Huang et al. (1998) extracts Intrinsic Mode Functions (IMFs) from a 1D signal by the sifting process which leaves the final residue as a constant or a monotone trend. An Intrinsic Mode Function (IMF) is originally defined that the number of its extrema equals the number of its zero-crossings and it has a zero local mean which represents a generally simple oscillatory mode as a counterpart to the simple harmonic function.
The first IMF is extracted from the given 1D signal g(t) by the following sifting process:
| Step 1: | Find all local extrema (minima and maxima) of g |
| Step 2: | Interpolate all local minima (resp. maxima) by the cubic spline to obtain the lower envelope LE (resp. the upper envelope UE) |
| Step 3: | Get the local mean m of the lower envelope LE and the upper envelope UE, i.e., m = (LE+UE)/2 |
| Step 4: | Compute h = g-m. If h satisfies the stopping criterion of the sifting process, then h is defined as an IMF, otherwise set g = h and repeat the process from step 1. The first IMF is denoted as f1 and is specified as the first residue |
The next IMF f2 is then extracted by applying the above sifting process (step 1-4) to the first residue r1. This process is repeated until the residue is a constant or a monotone trend. When this process is complete, the original signal g is decomposed as follows:
| (1) |
where, K is the number of IMFs f1, f2,..., fK and rK is the final residue. The leading IMFs contain high local spatial scales and/or variations, the trailing IMFs give coarse local spatial scales.
NEW BI-DIMENSIONAL EMD ALGORITHM
The essence of this new bi-dimensional EMD algorithm is to transform the given 2D data (e.g., an image) into a function defined on a triangular mesh and to use a novel interpolation method based on bi-Laplace operator to estimate envelopes which can accelerate the envelope computation of the traditional bi-dimensional EMD method based on thin plate spline interpolation (Linderhed, 2005). It can be summarized as follows:
| Step 1: | Transform the given 2D data into a function g defined on a triangular mesh M = (V, F) |
| Step 2: | Find all local extrema (minima and maxima) of the function g defined on the triangular mesh M = (V, F) |
| Step 3: | Interpolate all local minima (resp. maxima) by the method based on bi-Laplace operator to obtain the lower envelope LE (resp. the upper envelope UE) and get the local mean m of the lower envelope LE and the upper envelope UE, i.e., m = (LE+UE)/2 |
| Step 4: | Compute h = g-m. If h satisfies the stopping criterion of the sifting process, then h is defined as a BIMF (bi-dimensional intrinsic mode function), otherwise set g = h and repeat the process from step 2. The first BIMF is denoted as f1 and r1 = g-f1 is specified as the first residue |
The next BIMF f2 is then extracted by applying the above sifting process (step 2-4) to the first residue r1. This process is repeated until the residue is a constant or a monotone trend, or a specified number BIMFs have been extracted.
In the following, the key technique of the new bi-dimensional EMD algorithm will be introduced including transforming the given 2D data into a function defined on a triangular mesh, extremum definition, interpolation algorithm based on bi-Laplace operator and stop criteria.
Transform the given 2D data into a function defined on a triangular mesh: Without loss of generality, an MxN image g(m, n) can be used as the given 2D data and assume that the values g(m, n) correspond to an approximation of a continuous function defined on [0, 1]x[0, 1], at the set of points:
A triangular mesh M = (V, F) can be generated with the point set V and image boundary edges by the constrained Delaunay triangulation algorithm (Fig. 1a). V denotes the set of vertices of the mesh, F = {(i, j, k)|vi, vj are vk the three vertices of a triangle in M} denotes the set of faces.
Then the MxN image g(m, n) [m, n]ε{0,..., M-1}x{0,..., N-1} can be regarded as a function defined on the mesh M = (V, F). It is discretized as a piece-wise linear function, which is defined by linearly interpolating the values of g at the vertices, i.e., g = (g(v1), g(v2),..., g(vMxN)) using the barycentric coordinates (Fig. 1d-e).
Extremum definition: g(vi) is defined as a local minimum or a local maximum of the function g living on a triangular mesh, if it is lower or larger than all values of g at the 1-ring neighbors of the vertex vi.
| Fig. 1(a-e): | Transform 2D images into functions defined on triangular meshes. (a) An example for generating the triangular mesh of a 25x25 image, (b) Lena, (c) Camera man, (d) Visualization of the function defined on a triangular mesh for Lena and (e) Visualization of the function defined on a triangular mesh for camera man |
The set of vertex indices of the 1-ring neighbors of vi is N(i) = {j|vi, vj} is an edge of the mesh}.
Based on this definition, the indices of the local minima and maxima of g are represented as:
and:
Interpolation method based on bi-laplace operator: Thin plate spline is a natural extension of the cubic spline to 2D which is very important for data interpolation and approximation in Euclidean space. It can be characterized as the minimizer of a thin plate energy which refers to a physical analogy involving the bending of a thin sheet of metal (Bookstein, 1989).
Specifically, the thin plate energy of the function f defined on Ω⊂R2 is defined as:
| (2) |
The corresponding Euler-Lagrange equations is:
Δ2f = 0 |
(3) |
where, Δ2 is bi-Laplace operator. Therefore bi-Laplace operator minimizes the thin plate energy. Consequently, it can be used to generate an interpolation function subjecting to data interpolation constraints.
Consider the discretization of the 2D domain Ω as the triangular mesh M = (V, F), Laplace operator Δ can be discretized as a sparse matrix L with elements (Meyer et al., 2002):
| (4) |
where, |V| is the number of vertices of the mesh, N(i) is the set of vertex indices of the 1-ring neighbors of vi:
αij and βij are the angles opposite of the mesh edge (vi, vj) and Ai is the Voronoi area of vertex vi (Fig. 2).
| Fig. 2: | Discretization of Laplace operator |
Then the following optimization scheme based on bi-Laplace operator can be solved to interpolate all data constraints which is similar to the approach of Wang et al. (2012) and is an extension of the cubic spline to 2D Euclidean space:
| (5) |
where, L is the discretization Laplacian matrix 4, f = (f(v1), f(v2),..., f(v|V|)) is the unknown interpolation vector and C is the interpolated anchor set.
The optimization scheme 5 can lead to a sparse linear system for values of f at the vertices of the triangular mesh. Figure 3 gives the results for interpolating all local minima and maxima of Lena image by the method based on bi-Laplace operator to obtain the lower envelope and upper envelope.
| Fig. 3(a-d): | Envelope generation by the interpolation method based on bi-Laplace operator (a) Local minima, (b) Lower envelope interpolating all local minima, (c) Local maxima and (d) Upper envelope interpolating all local maxima |
Stop criteria: Just like the original 1D EMD, the stop criterion of the sifting process is controlled by limiting the size of the standard deviation SD which is computed from the values of the two consecutive sifting results hj-1 and hj at all vertices:
| (6) |
where, |V| denotes the number of the vertices.
The sifting process is stopped if SD falls below a threshold SDT. The typical value of SDT is set between 0.1 and 0.3 as the original EMD.
EXPERIMENT RESULT ANALYSIS AND DISCUSSION
The above bi-dimensional EMD algorithm was implemented on a computer with Inter Core(TM)2 Duo CPU T6660 at 2.20, 2.19 GHz and with 2.96 GB RAM. The backslash operator of MATLAB 2010 (Yan, 2011) is used to solve a sparse linear system. Although, the proposed method is capable of decomposing images of fine resolution fast, the maximum image size is limited to 256H256 pixel in order to compare it to the traditional EMD method based on thin plate spline interpolation (Linderhed, 2005).
An orthogonality index, denoted as OI, has been proposed for IMFs by Huang et al. (1998) to test the decomposition quality of an EMD algorithm which may be extended to the case of BIMFs for an MxN image g(m, n) as follows:
| (7) |
where, fi(x, y), i = 1,...K are BIMFs and fK+1(x, y) is the Kth residue. A low value of OI indicates a good decomposition in terms of local orthogonality among the BEMFs.
Figure 4 and 5 give the comparison results for Lena image and Camera Man image between the presented method and the traditional method based on thin plate spline interpolation (Linderhed, 2005), respectively. The values of BIMFs may be negative, so the values of all BIMFs are transformed into [0, 255] to visualize BIMFS. According to the comparison results (Fig. 4 and 5) and orthogonality statistics (Table 1) under the same conditions (SDT = 0.1, K = 3), the proposed method can generate comparable results (even better results) than the tradition algorithm based on thin plate spline interpolation. In addition, the presented method can accelerate the decomposition process greatly (Table 1).
| Table 1: | EMD comparison in orthogonality index (OI) and runtime between
the method of Linderhed (2005) and the presented method |
| Fig. 4(a-h): | EMD comparison for Lena between the method of Linderhed (2005) (a-d) and the presented method (e-h). (a)(e): BIMF1, (b)(f): BIMF2, (c)(g): BIMF3 and (d)(h): Third residue |
| Fig. 5(a-h): | EMD comparison for camera man between the method of Linderhed (2005) (a-d) and the presented method (e-h). (a)(e): BIMF1, (b)(f): BIMF2, (c)(g): BIMF3 and (d)(h): Third residue |
The main reason is that a sparse linear system is only solved when generating an envelope in the sifting process while the tradition algorithm based on thin plate spline interpolation (Linderhed, 2005) needs to solve a full linear system with the same number unknowns as data points.
Application in image enhancement: Image enhancement is one of major research areas in digital image processing, whose purpose is to bring out details that are hidden in an image, or to increase the contrast in a low contrast image. It is very beneficial to further image applications, such as image segmentation, image analysis and understanding, image classification and so on (Saradhadevi and Sundaram, 2010).
Bi-dimensional empirical mode decomposition can adaptively decompose an image g into several Bi-dimensional Intrinsic Mode Functions (BIMFs) fi, i = 1,...K with different scale details and a residue rK. In fact, the leading BIMFs contain the fine-scale details and the trailing BIMFs represent the smoothed features while the residue can be regarded as a monotone trend (i.e., the lightness of the image). Consequently, image details can be enhanced by modifying IMFs and image lightness can be improved by adjusting the residue. The specified adjustment formula based on the proposed bi-dimensional EMD algorithm is as follows:
| (8) |
where, Ti, i = 1,...K+1 are the transformation functions for different BIMFs and the residue, g’ represents an enhancement result for the input image. Generally, K = 3 and Ti, i = 1,..., K+1 selected as linear functions can generate good results (Fig. 6e), e. g.:
Ti = 1.8-0.2x(i-l), i = 1,
2, 3, 4 |
(9) |
According to the theory of Jobson et al. (2002), the overall lightness of an image can be measured by the image mean which is also the ensemble measure for regional lightness. The overall contrast can be measured by taking the mean of regional standard deviation which provides a gross measure of the regional contrast variations. The following image mean and variance changes can be used as a quantitative evaluation of an image enhancement algorithm:
| (10) |
| (11) |
where, g’ represents an enhancement result for an input MxN image g, Var(g’ (m, n)) and Var(g(m, n)) are regional standard deviations of g and g’ at the pixel (m, n) and Mean(.) represents the mean of an image.
| Fig. 6(a-e): | Image enhancement comparison for camera man (a) Original image with low contrast, (b) Histeq method, (c) Adapthisteq method, (d) Imadjust method and (e) The presented method with transformation functions 9 |
The above image enhancement method was compared to some classical methods, such as histogram equalization (Histeq), adaptive histogram equalization (Adapthisteq) and image contrast adjustment (Imadjust) (Chang and Wu, 1998; Gonzalez and Woods, 2006):
| • | Histeq performs histogram equalization which operates on the whole image in order that the histogram of the output image approximately matches the uniform distribution histogram |
| • | Adapthisteq performs contrast-limited adaptive histogram equalization which increases the contrast of the input image by operating on small data regions. The contrast of each data region is enhanced so that the histogram of each output region approximately matches the uniform distribution histogram |
| • | Imadjust improves the contrast of the image by mapping the values of the input intensity image to new values in order that a specified percent of the data (e.g., 1%) is saturated at low and high intensities of the input data |
| Table 2: | Image enhancement comparison in changes of contrast (C) and lightness (L) |
According to the comparison for camera man between the presented image enhancement method and these classical methods implemented by image processing toolbox of MATLAB 2010 (Yan, 2011), the presented image enhancement method can control better in contrast and lightness in image enhancement (Fig. 6) and obtain a better result about quantitative evaluation (Table 2).
CONCLUSION
In this study, a novel bi-dimensional EMD algorithm is presented for 2D data (e.g., images). The essence of this approach is based on transforming the given 2D data into a function defined on a triangular mesh and using an interpolation method based on bi-Laplace operator to generate envelopes. It can be implemented faster than the traditional bi-dimensional EMD method based on thin plate spline interpolation. Experiment results show that this method is robust and efficient and can be applied in image enhancement. The application examples also prove that the image enhancement algorithm based on the presented bi-dimensional EMD is efficient and better than some classical algorithms in contrast and lightness controlling.
In the future, the proposed bi-dimensional EMD algorithm can be used in other areas such as image denoising, image analysis, pattern recognition and so on.
ACKNOWLEDGMENTS
This study is supported by National Natural Science Foundation of China (No. 61202261, No. 61173102), Scientific and Technological Development Plan of Jilin Province (No. 20130522113JH) and opening Foundation of Key Laboratory of Symbolic Computation and Knowledge Engineering of Ministry of Education of China (No. 93K172012K02).