Yan-Li Liu et al.:Pores-Preserving Face Cleaning 559 transformation.The Fourier transform is designed to simplicity,we call it directional BEMD.The methods of work with liner and stationary signals.The wavelet the second category are genuine sense 2D EMD which transform,on the other hand,is well-suited to handle use radial basis functions(RBF)or cubic spline to in- non-stationary data,but is poor at processing non- terpolate envelopeslis-201.Although they have been linear data.However,few energy-frequency data are successfully used in some cases,these algorithms suffer truly linear and stationary.To address this problem from the problem of overshoot or undershoot because the Hilbert-Huang transform (HHT)has been recently of their interpolation scheme.In this paper,we pro- developedl4l.The HHT comprises two steps:EMD pose two solutions,namely limiting minimum frequency generates a finite number of intrinsic mode functions and adaptive local mean,to improve the performance (IMF),and then the Hilbert transform is applied to of EMD. IMF.The combination of IMF and its Hilbert trans- form forms an analytic signal,which can be used to 2.21 Background of EMD generate a "time-frequency-energy"representation of The principle of EMD is to identify the intrinsic os- the data. cillatory modes embedded in the signals based on their In this paper,we choose EMD and HHT rather than characteristic time scales,and then to decompose the wavelet transform due to the following considerations. signals into a sum of IMFs.In order to permit phys- First,imperfect facial images are in general nonlinear ically meaningful instantaneous frequency to be de- and nonstationary in nature,EMD is suitable to pro- fined over it,an IMF must satisfy the following two cess these data.Second,based on the local time scale conditionsl4l:first,the numbers of extrema and zero- of the data,EMD decomposes a signal into IMFs adap- crossings must differ by at most one;second,it is sym- tively without using a prior basis which do not neces- metric with respect to the local mean of the data.Ide- sarily match the varying nature of the signals.With- ally,the requirement of the second condition should be out the need of pre-specifying a decomposition level as "the local mean of the data is zero".Since it is difficult requested by wavelet decomposition to extract out all to define a local averaging time scale in general case,in the image's oscillatory modes,EMD allows us to use practice,the local mean of the data is usually replaced all the oscillatory modes to further define imperfect de- by the mean of the envelopes defined by interpolating gree.Thus,by exploiting EMD,we can obtain adap- the local maxima and local minima. tive multi-resolution representation of image.Due to EMD is an iterative process.To illustrate its con- unpredictable appearances of various kinds of distrac- cept,here we describe the algorithm of EMD in 1D tions,such an adaptive multi-resolution representation case.Given a signal x(t),the algorithm of EMD can be is favorable.Third,the IMF is generally in good agree- summarized as follows. ment with intuitive and physical signal interpretations, 1)Identify all the local extrema of x(t) and has well-defined instantaneous frequency,allowing 2)Interpolate between local minima (resp.local us to define local energy to capture the characteris- maxima),ending up with the envelope emin(t)(resp. tics of distractions and pores.On the contrary.the emax(t)). local energy extracted by Fourier or wavelet3-15]suf 3)Compute the mean m(t)=(emin(t)+emax(t))/2. fers from the problem of energy leakage,due to the 4)Extract the detail d(t)=x(t)-m(t).Check prop- infinite or limited finite length of the basis.The energy erty of d(t): leakage makes the quantitative definition of the energy- if it meets the above two conditions,an IMF frequency-space distribution of image difficult is derived,and replace x(t)with the residual r(t)= Due to the nature of image data,when process- x(t)-d(t): ing images with idea of EMD,bidimensional empir- if it is not an IMF,replace r(t)with d(t). ical mode decomposition (BEMD)is necessary.Ac- 5)Repeat Steps 1)~4)until the residual satisfies cordingly,to perform 2D space-frequency analysis of some stopping criterion. IMFs,Riesz transform,an isotropic generalization of Finally,x(t)is represented as the sum of a finite the Hilbert transform to multiple dimensions is also number of IMFs (both amplitude and frequency mod- needed[16).For BEMD,existing methods can be di- ulated)and a residual trend.The process of iteratively vided into two categories.The first category of the extracting an IMF is also called sifting process.One can approaches divide an image into one-dimensional data. observe from the above formulation that there is no pre- Then the 1D EMD is applied to a limited number of ori- fixed basis involving in the sifting process,instead,the entations:two (horizontal and vertical)or morel171.For decomposition proceeds depending on the data itself.Yan-Li Liu et al.: Pores-Preserving Face Cleaning 559 transformation. The Fourier transform is designed to work with liner and stationary signals. The wavelet transform, on the other hand, is well-suited to handle non-stationary data, but is poor at processing non￾linear data. However, few energy-frequency data are truly linear and stationary. To address this problem, the Hilbert-Huang transform (HHT) has been recently developed[4]. The HHT comprises two steps: EMD generates a finite number of intrinsic mode functions (IMF), and then the Hilbert transform is applied to IMF. The combination of IMF and its Hilbert trans￾form forms an analytic signal, which can be used to generate a “time-frequency-energy” representation of the data. In this paper, we choose EMD and HHT rather than wavelet transform due to the following considerations. First, imperfect facial images are in general nonlinear and nonstationary in nature, EMD is suitable to pro￾cess these data. Second, based on the local time scale of the data, EMD decomposes a signal into IMFs adap￾tively without using a prior basis which do not neces￾sarily match the varying nature of the signals. With￾out the need of pre-specifying a decomposition level as requested by wavelet decomposition to extract out all the image’s oscillatory modes, EMD allows us to use all the oscillatory modes to further define imperfect de￾gree. Thus, by exploiting EMD, we can obtain adap￾tive multi-resolution representation of image. Due to unpredictable appearances of various kinds of distrac￾tions, such an adaptive multi-resolution representation is favorable. Third, the IMF is generally in good agree￾ment with intuitive and physical signal interpretations, and has well-defined instantaneous frequency, allowing us to define local energy to capture the characteris￾tics of distractions and pores. On the contrary, the local energy extracted by Fourier or wavelet[13−15] suf￾fers from the problem of energy leakage, due to the infinite or limited finite length of the basis. The energy leakage makes the quantitative definition of the energy￾frequency-space distribution of image difficult. Due to the nature of image data, when process￾ing images with idea of EMD, bidimensional empir￾ical mode decomposition (BEMD) is necessary. Ac￾cordingly, to perform 2D space-frequency analysis of IMFs, Riesz transform, an isotropic generalization of the Hilbert transform to multiple dimensions is also needed[16]. For BEMD, existing methods can be di￾vided into two categories. The first category of the approaches divide an image into one-dimensional data. Then the 1D EMD is applied to a limited number of ori￾entations: two (horizontal and vertical) or more[17]. For simplicity, we call it directional BEMD. The methods of the second category are genuine sense 2D EMD which use radial basis functions (RBF) or cubic spline to in￾terpolate envelopes[18−20]. Although they have been successfully used in some cases, these algorithms suffer from the problem of overshoot or undershoot because of their interpolation scheme. In this paper, we pro￾pose two solutions, namely limiting minimum frequency and adaptive local mean, to improve the performance of EMD. 2.2 Background of EMD The principle of EMD is to identify the intrinsic os￾cillatory modes embedded in the signals based on their characteristic time scales, and then to decompose the signals into a sum of IMFs. In order to permit phys￾ically meaningful instantaneous frequency to be de- fined over it, an IMF must satisfy the following two conditions[4]: first, the numbers of extrema and zero￾crossings must differ by at most one; second, it is sym￾metric with respect to the local mean of the data. Ide￾ally, the requirement of the second condition should be “the local mean of the data is zero”. Since it is difficult to define a local averaging time scale in general case, in practice, the local mean of the data is usually replaced by the mean of the envelopes defined by interpolating the local maxima and local minima. EMD is an iterative process. To illustrate its con￾cept, here we describe the algorithm of EMD in 1D case. Given a signal x(t), the algorithm of EMD can be summarized as follows. 1) Identify all the local extrema of x(t). 2) Interpolate between local minima (resp. local maxima), ending up with the envelope emin(t) (resp. emax(t)). 3) Compute the mean m(t) = (emin(t) + emax(t))/2. 4) Extract the detail d(t) = x(t)−m(t). Check prop￾erty of d(t): • if it meets the above two conditions, an IMF is derived, and replace x(t) with the residual r(t) = x(t) − d(t); • if it is not an IMF, replace x(t) with d(t). 5) Repeat Steps 1)∼4) until the residual satisfies some stopping criterion. Finally, x(t) is represented as the sum of a finite number of IMFs (both amplitude and frequency mod￾ulated) and a residual trend. The process of iteratively extracting an IMF is also called sifting process. One can observe from the above formulation that there is no pre- fixed basis involving in the sifting process, instead, the decomposition proceeds depending on the data itself