INTRODUCTION
EMD (Empirical mode decomposition) is developed by Huang
et al. (1998) and Huang and Wu (2008) which
can adaptively decompose a 1dimensional (1D) nonlinear and nonstationary 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 2dimensional (2D) data (e.g., images) by
EMD technique, i.e., bidimensional EMD (BEMD). Currently, several methods about
bidimensional 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 bidimensional
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 finiteelement 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 bidimensional
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 bidimensional EMD
algorithm using thin plate spline interpolation based on biLaplace 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 zerocrossings 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 = gm. 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 f_{1} and is specified
as the first residue 
The next IMF f_{2} is then extracted by applying the above sifting
process (step 14) to the first residue r_{1}. 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:
where, K is the number of IMFs f_{1}, f_{2},..., f_{K}
and r_{K} is the final residue. The leading IMFs contain high local
spatial scales and/or variations, the trailing IMFs give coarse local spatial
scales.
NEW BIDIMENSIONAL EMD ALGORITHM
The essence of this new bidimensional 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 biLaplace operator to estimate envelopes
which can accelerate the envelope computation of the traditional bidimensional
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 biLaplace
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 = gm. If h satisfies the stopping criterion of the sifting
process, then h is defined as a BIMF (bidimensional intrinsic mode function),
otherwise set g = h and repeat the process from step 2. The first BIMF is
denoted as f_{1} and r_{1} = gf_{1} is specified
as the first residue 
The next BIMF f_{2} is then extracted by applying the above sifting
process (step 24) to the first residue r_{1}. 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 bidimensional 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 biLaplace 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)v_{i},
v_{j} are v_{k} the three vertices of a triangle in M} denotes
the set of faces.
Then the MxN image g(m, n) [m, n]ε{0,..., M1}x{0,..., N1} can
be regarded as a function defined on the mesh M = (V, F). It is discretized
as a piecewise linear function, which is defined by linearly interpolating
the values of g at the vertices, i.e., g = (g(v_{1}), g(v_{2}),...,
g(v_{MxN})) using the barycentric coordinates (Fig. 1de).
Extremum definition: g(v_{i}) 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 1ring neighbors of the vertex v_{i}.

Fig. 1(ae): 
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 1ring neighbors of v_{i} is N(i)
= {jv_{i}, v_{j}} 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 bilaplace 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 Ω⊂R^{2}
is defined as:
The corresponding EulerLagrange equations is:
where, Δ^{2} is biLaplace operator. Therefore biLaplace 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):
where, V is the number of vertices of the mesh, N(i) is the set of vertex
indices of the 1ring neighbors of v_{i}:
α_{ij} and β_{ij} are the angles opposite of the
mesh edge (v_{i}, v_{j}) and A_{i} is the Voronoi area
of vertex v_{i} (Fig. 2).

Fig. 2: 
Discretization of Laplace operator 
Then the following optimization scheme based on biLaplace 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:
where, L is the discretization Laplacian matrix 4, f = (f(v_{1}), f(v_{2}),...,
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 biLaplace operator to obtain the lower envelope and upper envelope.

Fig. 3(ad): 
Envelope generation by the interpolation method based on biLaplace
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
h_{j1} and h_{j} at all vertices:
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 bidimensional 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:
where, f_{i}(x, y), i = 1,...K are BIMFs and f_{K+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(ah): 
EMD comparison for Lena between the method of Linderhed
(2005) (ad) and the presented method (eh). (a)(e): BIMF1, (b)(f):
BIMF2, (c)(g): BIMF3 and (d)(h): Third residue 

Fig. 5(ah): 
EMD comparison for camera man between the method of Linderhed
(2005) (ad) and the presented method (eh). (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).
Bidimensional empirical mode decomposition can adaptively decompose an image
g into several Bidimensional Intrinsic Mode Functions (BIMFs) f_{i},
i = 1,...K with different scale details and a residue r_{K}. In fact,
the leading BIMFs contain the finescale 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 bidimensional EMD algorithm
is as follows:
where, T_{i}, 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 T_{i}, i = 1,..., K+1 selected as linear
functions can generate good results (Fig. 6e), e. g.:
T_{i} = 1.80.2x(il), 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:
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(ae): 
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 contrastlimited 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 bidimensional 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 biLaplace operator to generate envelopes. It can be implemented
faster than the traditional bidimensional 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 bidimensional EMD
is efficient and better than some classical algorithms in contrast and lightness
controlling.
In the future, the proposed bidimensional 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).