560 J.Comput.Sci.Technol.,May 2009,Vol.24.No.3 3 Improved EMD Sk=2k+1+1. .Initialize hk.o(p)=rk-1(p),11. Theoretically,there are two factors incurring the drawback of classic EMD mentioned in the above Detect all local extreme points of hk_1.Estab- lish the local maximum point set and local minimum section.First,since there is no bandwidth constraint during the decomposition process in classical EMD, point set respectively.For every pixel p of hk.1,if over-iterating on the whole signal for a better local ap- I(p)is strictly higher or lower than all its 8-connected neighbors,it is a local extreme point.Otherwise,it is proximation has the drawback of contaminating other parts of the signal,in particular,in uniformizing the not a local extreme point. Compute the geometric distance (Euclidean dis- amplitude and in over-decomposing it by spreading out its components over adjacent modes21.Second,clas- tance)between every two adjacent local maximum sical EMD uses the mean of envelopes interpolated by points row-wisely and then column-wisely in For the local maxima and minima respectively as a sub- example,suppose two adjacent local maximum points stitute for the local mean of data.Unfortunately,the in the same column are (1,y),(2,y)and 1 2, interpolation method,for example spline fitting,often then the distance S between them is z2-z1.If S>S, overshoots or undershoots the strongly nonlinear and then put the pixels(1+ix S,y),i=1,2,...,(S/S)] nonstationary data.Such overshoots and undershoots into The purpose of supplementing is to make the distance between every two adjacent local maxi- would be amplified through many times of spline fitting in the sifting process,and eventually result in diver- mum points not more than S,ensuring the distribution gence of IMF.Therefore,unless a better approximation of extreme points is not too sparse.Furthermore,to re- duce the boundary effects,we add extreme surrounding to local mean is conducted,EMD cannot be improved significantly. the original image boundaries by mirror symmetry.The For each of the above two problems,we propose a so- same operations are also performed on. For every pixel p,compute its local mean lution to get around it.The first solution is related with time-frequency uncertainty principle.As it is known, emean,1-1(p)using the adaptive local mean algorithm de- scribed below. uncertainty principle formulated as ·hk,l(p）=hk,l-1(p)-Cmean,.l-1(p): 7f×7t≥1/2 (1) If the stop criterion is met4,imf(p)=hk.(p), and stop the sifting process;else I=I+1. defines the relation between the width of frequency win- 3)rk(p)=rk-1(p)-imfk(p). dow Vf and that of the time window Vt.According 4)If the r(p)is a monotonic or constant signal, to (1),if we set Vt s a (a 0),then it must have terminate the decomposition process.Otherwise go to Vf >b(b=1/2a).That is,if we limit the maximal Step 2). temporal resolution,a minimal frequency resolution can The image I is finally expressed as the sum of finite be confined correspondingly.In the paper,we specify IMFs and a residue,as shown in the following equation: a maximum neighborhood in sifting process,as a re- sult,the minimal frequency of the resulting IMF is de- fined.Our solution to the second problem is:instead of I(p)= im时fk(p)+rK(p) (2) k=1 endeavoring to develop a fitting method which always cause overshoot or undershoot,we directly use the local The adaptive local mean algorithm is shown in Ta- mean of the image data to ascertain symmetry.To do ble 1:the function NumErtrema(p)returns the num- so,the only issue is to define a local space scale.Fortu- ber of local extreme points in p.r.Symmetry(p) nately,the above specified maximum neighborhood is tests whether the local extreme points in p.r are spa- a good candidate for defining such a local space scale. tially symmetric,and is defined as follows.The col- Taking computation cost and locality into account,an umn and row of p equally divide pr into four quad- adaptive local mean algorithm is developed. rants.Symmetry(p)is true if the following two condi- Integrating the above two strategies into one frame- tions are both satisfied:first,in each of the four quad- work.for a given image I with size of M x N in pixels, rants one local extreme point can be found;second, our ALNEMD algorithm is formulated as the following among these four local extreme points,two of them are steps. local maximum points and the other two are local min- 1)Initialize ro(p)=I(p). imum points.Otherwise Symmetry(p)is false.The 2)Extract k-th IMF imf&. threshold of ur is used to guarantee that the number ·Set the size of maximum neighborhood Sx×Sk, of local extreme points is not too few and typically set560 J. Comput. Sci. & Technol., May 2009, Vol.24, No.3 3 Improved EMD Theoretically, there are two factors incurring the drawback of classic EMD mentioned in the above section. First, since there is no bandwidth constraint during the decomposition process in classical EMD, over-iterating on the whole signal for a better local ap￾proximation has the drawback of contaminating other parts of the signal, in particular, in uniformizing the amplitude and in over-decomposing it by spreading out its components over adjacent modes[21]. Second, clas￾sical EMD uses the mean of envelopes interpolated by the local maxima and minima respectively as a sub￾stitute for the local mean of data. Unfortunately, the interpolation method, for example spline fitting, often overshoots or undershoots the strongly nonlinear and nonstationary data. Such overshoots and undershoots would be amplified through many times of spline fitting in the sifting process, and eventually result in diver￾gence of IMF. Therefore, unless a better approximation to local mean is conducted, EMD cannot be improved significantly. For each of the above two problems, we propose a so￾lution to get around it. The first solution is related with time-frequency uncertainty principle. As it is known, uncertainty principle formulated as ∇f × ∇t > 1/2 (1) defines the relation between the width of frequency win￾dow ∇f and that of the time window ∇t. According to (1), if we set ∇t 6 a (a > 0), then it must have ∇f > b (b = 1/2a). That is, if we limit the maximal temporal resolution, a minimal frequency resolution can be confined correspondingly. In the paper, we specify a maximum neighborhood in sifting process, as a re￾sult, the minimal frequency of the resulting IMF is de- fined. Our solution to the second problem is: instead of endeavoring to develop a fitting method which always cause overshoot or undershoot, we directly use the local mean of the image data to ascertain symmetry. To do so, the only issue is to define a local space scale. Fortu￾nately, the above specified maximum neighborhood is a good candidate for defining such a local space scale. Taking computation cost and locality into account, an adaptive local mean algorithm is developed. Integrating the above two strategies into one frame￾work, for a given image I with size of M × N in pixels, our ALNEMD algorithm is formulated as the following steps. 1) Initialize r0(p) = I(p). 2) Extract k-th IMF imf k . • Set the size of maximum neighborhood Sk × Sk, Sk = 2k+1 + 1. • Initialize hk,0(p) = rk−1(p), l = 1. • Detect all local extreme points of hk,l−1. Estab￾lish the local maximum point set A and local minimum point set B respectively. For every pixel p of hk,l−1, if I(p) is strictly higher or lower than all its 8-connected neighbors, it is a local extreme point. Otherwise, it is not a local extreme point. • Compute the geometric distance (Euclidean dis￾tance) between every two adjacent local maximum points row-wisely and then column-wisely in A . For example, suppose two adjacent local maximum points in the same column are (x1, y), (x2, y) and x1 < x2, then the distance S between them is x2 −x1. If S > Sl , then put the pixels (x1+i×Sl , y), i = 1, 2, . . . , b(S/Sl)c into A . The purpose of supplementing A is to make the distance between every two adjacent local maxi￾mum points not more than Sl , ensuring the distribution of extreme points is not too sparse. Furthermore, to re￾duce the boundary effects, we add extreme surrounding the original image boundaries by mirror symmetry. The same operations are also performed on B. • For every pixel p, compute its local mean emean,l−1(p) using the adaptive local mean algorithm de￾scribed below. • hk,l(p) = hk,l−1(p) − emean,l−1(p). • If the stop criterion is met[4] , imf k (p) = hk,l(p), and stop the sifting process; else l = l + 1. 3) rk(p) = rk−1(p) − imf k (p). 4) If the rk(p) is a monotonic or constant signal, terminate the decomposition process. Otherwise go to Step 2). The image I is finally expressed as the sum of finite IMFs and a residue, as shown in the following equation: I(p) = X K k=1 imf k (p) + rK(p). (2) The adaptive local mean algorithm is shown in Ta￾ble 1: the function NumExtrema(p) returns the num￾ber of local extreme points in Ωp,T . Symmetry(p) tests whether the local extreme points in Ωp,T are spa￾tially symmetric, and is defined as follows. The col￾umn and row of p equally divide Ωp,T into four quad￾rants. Symmetry(p) is true if the following two condi￾tions are both satisfied: first, in each of the four quad￾rants one local extreme point can be found; second, among these four local extreme points, two of them are local maximum points and the other two are local min￾imum points. Otherwise Symmetry(p) is false. The threshold of µT is used to guarantee that the number of local extreme points is not too few and typically set