Proceedings of the 6th World Congress on Intelligent Control and Automation, June 21 - 23, 2006, Dalian, China

Image =ooming Method Using 2D EMD Technique Yan Tian, Ying Huang, Yongjiang Li Department of Electronics & Information Engineering Huazhong University of Science and Technology Wuhan, 430074, China {tianyan2000, huangying581, liyongjiang119}@126.com Abstract - Image zooming is a prime and important problem in image processing. An image decomposition method called empirical mode decomposition (EMD) is first addressed in this paper, and then the problems such as searching local extrema points and boundary extension are discussed in detail. On this basis, a new image zooming method using 2 dimensional (2D) EMD technique is proposed. The main procedures are designed as the following: first, by 2D EMD method, an image is decomposed into several layers, where each of them contains different details or possesses different roughness; second, each layer is zoomed using a combination method with different parameters that are proportional to roughness of that layer, at last, the final zooming result is reconstructed from the combined interpolated layers. Experimental results demonstrate the effectiveness of the proposed method. Index Terms - Image zooming, 2 dimensional EMD, Combination interpolation

I. INTRODUCTION In order to acquire more details from an image to precisely identify and translate it, certain resolution must be assured for that image. There are two measures to increase the resolution: one is to improve the quality of sensors, the other is to design certain effective interpolation algorithms in software point of view. The former is bounded by hardware technique and economic factors generally. Thus, a reasonable zooming algorithm to increase the resolution is extremely necessary. Image zooming has gained important applications in various circumstances, such as improving remote sensing image spatial resolution, CCD picture zooming, high resolution medical imaging like CT, MRI, and ultrasonic imaging devices and systems, and so on. Aiming at various application fields, people have put forward several interpolation methods. Using an example sensor model, Price and Hayes gave out an optimized pre-filter for image interpolation [1]. In 2001, Wang and Ward made an interpolation method based on edge direction and enhancement [2]. Bai and Wang invented a dynamic online fuzzy error interpolation technique for robot compensation orientation [3]. Using the transmit equation, Battiato and Gallo provided an interpolation method based on anisotropy diffusing [4]. In the point of the local activity of images, Hadhoud’s researched a warped adaptive distance method for image interpolation [5]. However, in practice, the preciseness and complexity of algorithm need to be considered together.

1-4244-0332-4/06/$20.00 ©2006 IEEE

Therefore, several classical methods, such as bilinear interpolation, bi-cubic interpolation, and spline interpolation, are still widely used recently. However, none of these methods utilized the correlation between pixels, and thus sometimes the zooming effect is not satisfying. Empirical mode decomposition (EMD) is a data analysis method that is first introduced by Norden E. Huang in 1998 [6]. As a brand new time-frequency domain analytical technique, it is independent with Fourier methods, and is adaptive for analyzing nonlinear and non-stationary data sets. The concept of EMD is to decompose the signal into a set of functions defined by the signal itself, called Intrinsic Mode Functions (IMFs), which preserve the inherent properties of the original signal. As the increasing of the decomposition level, the complexity of the IMFs decreases, and so does the scale of the signal. Thus, different IMF reveals the state of the signal in varied scales as well as the internal relationship between the points over various ranges. This technique has been widely applied in 1D signal processing, while its application in image processing field is rather limited. This paper proposes a new image zooming method based on 2D EMD technique. After an image is decomposed by the EMD process, the decomposed layers are zoomed separately using a combinational method of bilinear and bi-cubic interpolation. The combination parameter is defined by the roughness of the corresponding layer. At last, the final result is reconstructed from the expanded image layers. This method distinguishes the shades among the layers by accounting for interrelationship of the pixels locally and globally. Therefore, it overcomes the shortcomings of traditional methods to some extend. II. 2D EMD It is known to all, wavelet transform is more adaptable for non-stationary signal analysis. Virtually, wavelet is a windowed Fourier transform, and finite length of wavelet basis may cause energy linkage. Moreover, once the wavelet base and the decomposition scale is given, the resultant signal is restricted to a fixed frequency band, which is solely dependent on sampling frequency and has no relationship with the signal itself. In this sense, wavelet is not adaptive. Compared with Wavelet transform, EMD shows a superior performance on selectivity and precision of data analysis. It is a powerful tool for adaptive multi-scale analysis of short time nonlinear and non-stationary signals. These interesting aspects

10036

of the EMD motivate the extension of this method to image processing. In this section, we first introduce the EMD process for images and solve some key problems including extrema location and boundary extension, and then utilize it in image zooming. A. General Decomposition procedure As mentioned above, EMD method decomposes a signal into a set of components called IMF, composing of a single oscillation mode, and a residue term. In 1D case, IMF components must satisfy two conditions [6]: 1) in the whole data set, the number of extrema and the number of zero crossings must either equal or differ at most by one; 2) at any point, the mean value of the envelopes defined by the local maxima and minima is zero. IMF is usually obtained through an iterative way called “sifting process”. Before we extract an IMF from the signal, a sifting process may take several rounds to refine it. A round of sifting process is described as: Step 1: find all local maximum and minimum of the current mode signal hl ( k −1) ( x, y ) , where hl ( k −1) ( x, y ) denote for the (k – 1) st iteration process sifting for the l-th IMF; and initially, h00 ( x, y ) is defined as the original image f ( x, y ) . Step 2: use proper interpolation method to build the upper envelop emax ( x, y ) and lower envelop emin ( x, y ) which pass through all local maximal and minimal points respectively. Step 3: subtract the hlk ( x, y ) with the mean envelop mlk ( x, y ) =

emax ( x, y ) + emin ( x, y ) , and the result is denoted as 2

hlk ( x, y ) = hl ( k −1) ( x, y ) − mlk ( x, y ) . If it is not an IMF signal, then the algorithm go to Step 1, until the criteria are satisfied. In the former statement, the two basic criteria seem too strict to extract an IMF in 2-dimensional case. Practically, we select an approximate IMF with close number of local maxima and minima, and the mean value of the envelopes defined by the local maxima and minima is close to zero. Hence, the criterion to stop the sifting process becomes: mlk ( x, y ) < ε , ∀( x, y ) . (1) Actually, a more reasonable method is to use two following conditions: mlk ( x, y ) < ε 1 , ∀( x, y ) (2) P{ mlk ( x, y ) > ε 2 } < θ .

(3)

Where ε1 and ε 2 are a larger and smaller threshold respectively. When the sifting process is ended, the updated mode function hlk ( x, y ) is an IMF component of the original signal. And it is denoted as: f (l ) ( x, y ) = hlk ( x, y ) . (4) Then, the residue r (l ) ( x, y ) is

r (l ) ( x, y ) = hl 0 ( x, y ) − hlk ( x, y ) . (5) If the residue still has many local maxima and minima, then it calls for further decomposition taking on the residue. Then it is needed to start a new sifting process for l + 1st IMF, h(l +1)0 ( x, y ) = r ( l ) ( x, y ) . (6) Continue extracting other IMF components while the amount of local maxima and minima excesses the required number. When the procedure finishes, the reminding residue r ( n ) ( x, y ) is the reminder of the decomposition. Thereby, image f ( x, y ) is decomposed to n – 1 IMF and the residue term r ( n ) ( x, y ) , that is n −1

f ( x, y ) = ¦ f (l ) ( x, y ) + r ( n ) ( x, y ) .

(7)

l =1

What’s more, in the actual decomposition procedure, there still exist some problems needed to be resolved, such as how to interpolate the scattered data so as to ensure the sifting process be convergent; how to search all local maxima and minima; and how to treat with the boundary interference and predict the local extrema outside the boundary. In the following sections, we will discuss these problems respectively. B. Interpolation method 2D EMD technique highly depends on a reasonable interpolation method. The number of local extrema must decrease; otherwise, the decomposition procedure cannot be convergent. So, in implementation of 2D EMD, one of key problems is how to fit a surface to the 2D scattered data points representing the local extrema. Generally, the solution requires the interpolated surface be smooth enough and go through every data points, overshot shall be avoided and the second derivative must be continuous everywhere for the signal. By these means, few new extrema points will be generated and the number of the existed local extrema may decline. The literature [7] introduced two interpolation methods including triangle-based cubic spline interpolation and a radial basis function called thin-plate smoothing spline interpolation. The difference between them relies on the general smoothness and computational cost. When it is to fit a surface to a few data points, the latter method generates a smoother surface with no “crease effect”. While the surface produced by triangle-based method is just piece-wise continuous. However, thin-plate smoothing spline interpolation is really a time-consuming process. By this method, any point we find depends on all given points, but actually, some points may have no connection at all. Thus, the computation cost soars. On the contrary, given a large pool of the data set, the fitted surface can be smooth enough by triangle-based cubic spline interpolation. And since it only takes on the neighboring points while neglecting all points that are in far away from the center point, the computational cost is much lower. For these reasons, in this paper we employ triangle-based cubic spline interpolation in 2D EMD process. C. Extrema point locating process

10037

It is easy to locate all extrema of a one-dimensional signal by searching and comparing along one direction. However, in two-dimensional case, the same scheme may fail to work due to saddle points. An image is, in fact, a set of 2D discrete array, and thus the extrema points are those whose gray values are larger or smaller than that of the neighboring points’. Thereby, a searching algorithm is designed as following: Case 1: if the center point is larger than all neighbor points, then the center point is a local maximum; Case 2: if the center point is smaller than all neighbor points, and then the center point is a local minimum; Case 3: if the center point is larger than some neighbor points while smaller than other points at the same time, then the center point is neither a local maximum nor minimum; Case 4: if the center point is equal to the extrema points in the neighborhood, then combine them together to form a region and continue comparing the region and its neighborhood until it is in one of former 3 cases or the growing region reaches the border. Because this algorithm compares the center point with all its neighbor points rather than comparing along certain direction, it avoids taking the saddle points as local extrema, and thus locate all local extrema points accurately. D. Boundary extension To avoid boundary effect, the extrema points should be extended. For 1D a signal, the common used methods include symmetry and linear prediction. Zeng introduced a simple edge extrapolation method, whose basic thought is to do even extension and odd extension for data [8]. We employ this idea to 2D case. The process of extrapolation is shown in the figure below. An image is reflected along four edges and symmetrized about the four corner points as well. In this way the original image is enlarged by nine times, and all the boundaries are ensured to be outside the processing area.

(a) original (b) extended Fig. 1 Boundary extension for bi-dimensional image.

In practice, searching work for the local extreme points are referred to of the original image, and the extension is done just for new extreme point coordinates, as a result, the computation of the searching algorithm would not be increased. III. AN IMAGE ZOOMING METHOD BASED ON 2D EMD

An image zooming problem can be regarded as a filter problem. In design, we often try to make the frequency response of the filter to approximate the frequency of an image as closely as possible [9]. Bilinear and bi-cubic are two widely used image zooming methods. Bilinear method utilizes the 4 nearest neighbor pixels and a new pixel is created by averaging the neighbor points. Similarly, bi-cubic makes use of the 16 nearest neighbor points and it is a higher order interpolation method. If an image contains much lower frequency information, it is reasonable to apply the bilinear interpolation than the bi-cubic method. On the contrary, it is suitable to use bi-cubic. The bilinear and bi-cubic formulae are given as follows: f1 ( x, y ) = a0 + a1 x + a2 y + a3 xy , (8) f 2 ( x, y ) = b0 + b1 x + b2 y + b3 x 2 + b4 xy

. (9) +b5 y 2 + b6 x 3 + b7 x 2 y + b8 xy 2 + b9 y 3 Generally, an image may contain both low frequency and high frequency ingredients. Therefore, it is reasonable to generate a new method that combines the advantages of bilinear and bi-cubic. Based on above discussions, the new image zooming method is designed as the following: Step 1: decompose the original image f ( x, y ) into n − 1 layers f (i ) ( x, y ) , (i = 1, 2" n − 1) and the residue term r ( n ) ( x, y ) using EMD method; Step 2: use combinational method to interpolate in each layer, and the zoomed image layers can be expressed as: F (i ) ( x, y ) = (1 − ρ ) f1( i ) ( x, y ) + ρ f 2( i ) ( x, y ) , (10) R ( n ) ( x, y ) = (1 − ρ )r1( n) ( x, y ) + ρ r2( n ) ( x, y ) .

(11)

Where i = 1, 2" n − 1 , F ( x, y ) is the corresponding layer (i )

of i-th layer IMF component. R ( n ) ( x, y ) is the residue term of the zoomed image. f1(i ) ( x, y ) and f 2 (i ) ( x, y ) are the interpolated result of the i-th layer produced by the bilinear and bi-cubic method respectively. r1( n ) and r2( n ) are the interpolated residue terms produced by the bilinear and bi-cubic method . ρ is the combination parameter. Step 3: based on the interpolated results in each layer, apply the formula (7) to reconstruct the final zoomed image. The combination parameter determines the weight of bi-cubic method, from the above analysis, for a complex image, ρ need to be increased; on the contrary, for a plain image, ρ should decrease. Since an image is not always uniform, the parameter ρ should not be fixed for a whole image. In this sense, we decompose the image into different layers and apply a different parameter ρi in each layer. Here, ρi is closely related with the complexity of the corresponding layer. Because the variance is a measure of the complexity and roughness of an image, it is convenient to choose a ρi according to the weight of the variance the in that layer, that is,

10038

ρi =

σ i2

, i = 1, 2" n − 1 .

n −1

¦σ j =1

(12)

this sense, IMFs reveal the multi-scale features of the image and the local and global characteristics as well.

2 j

Where, n is the number of layer after EMD procedure. The n -th level residue is usually a monotonous plane. Thus, it is suitable to use bilinear interpolation method and need not to count its variance, i.e., ρn = 0 . Other layers are interpolated in the combinational way. IV. EXPERIMENTAL RESULT AND ANALYSIS In this section, we design two groups of experiments. The first is aiming at demonstrate the EMD process for an image, and the second is to check zooming effect of the proposed method. A. Decomposition results The standard image Cameraman and Terrain, shown in Fig. 2, are chosen as the test data, Here, we choose the thresholds mentioned in formula (2) and (3) as ε1 = 100 , ε 2 = 20 , θ = 0.2 . The decomposed layers are shown in Fig. 3 and Fig. 4.

B. Image zooming result In order to verify the effectiveness of the proposed method in image zooming, we first interpolated the obtained layers using our method, and then recover the zoomed image from them. Fig. 5 (a) and (b) show the comparing experiments for Cameraman that apply bi-cubic and bilinear directly. Fig. 5 (c) shows the result of our method. The combination parameters we choose by formula (12) are 0.3849, 0.3509, 0.1672, 0.0970, 0 (the last layer use bilinear method directly). Fig. 6 (a) and (b) show the comparing experiments for Terrain. Fig. 6 (c) shows the result of our method. The combination parameters we choose by formula (12) are 0.4469, 0.4291, 0.1240, 0 (the last layer use bilinear method directly).

(a) Bicubic method

(b) Bilinear method

(c) Our method

Fig. 5 The zooming results for Cameraman. (a)

(b)

Fig. 2 Experiment data: (a) Cameraman; (b) Terrain.

(a)

(b)

(c)

(d)

(e)

(a) Bicubic method

(b) Bilinear method

(c) Our method

Fig. 6 The zooming results for Terrain.

(f)

(g)

(h)

(i)

(j)

Fig. 3 Decomposition results for Cameraman: (a)~(d) are IMF1~IMF4, (e) is the residue term, (f)~(j) are the corresponding residue for the 1st to 5-th decomposition layer respectively.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

From Fig. 5 (a), (b) and (c), it is obvious that (c) is clearest in the highlighted area, while (a) and (b) are dim and fuzzy. The similar effect can be observed from Fig. 6. So compared with bi-cubic and bilinear method, our combinational method that based on EMD technique has superior performance in vision effect. To assess the virtue of our method objectively, we calculate the PSNR (peek noise-to-signal ratio) and RMSE (root of mean square error) for the above results. From Tab. I, it is easy to draw a conclusion that our method is better than bi-cubic and bilinear by means of PSNR and RMSE.

Fig. 4 Decomposition results for Terrain: (a)~(c) are IMF1~IMF3, (d) is the residue term, (e)~(h) are the corresponding residue for the 1st to 4-th decomposition layer respectively.

From Fig. 3 and Fig. 4, it is clear that the complexity of each IMF decreases, and the residue term is almost a plane. In

10039

TABLE I PSNR AND RMSE OF ZOOMING RESULTS Bicubic Bilinear Our method PSNR 50.7662 50.8386 50.8829 Cameraman RMSE 19.276 19.2063 19.1638 PSNR 67.6803 67.7852 67.8904 Terrain RMSE 7.3926 7.3539 7.3154

V. CONCLUSION This paper first discusses the 2D EMD process and gives practical solutions to the problem of searching local extrema and boundary extension. Then a novel image zooming method based on 2D EMD is presented. The new method applies a combinational interpolation scheme in every decomposed image layers. Comparing with bi-cubic and bilinear methods, the experimental results demonstrate that the proposed method is effective and promising. Considering the complexity of EMD algorithm, our future researches focus on exploring a fast algorithm for 2D EMD and its applications in other image processing fields. ACKNOWLEDGMENT The work described in this paper was supported by National Natural Science Foundation of China (60272099) and Natural Science Foundation of Hubei (2004ABA051). The authors would like to thank Prof. Jian Liu for his valuable suggestions on this paper. REFERENCES [1]

[2] [3]

[4]

[5]

[6]

[7]

[8]

[9]

Price, J.R., Heyes, M.H., Optimal Pre-filtering for Improved Image Interpolation. Conference Record of the Thirty-Second Asilomar Conference on Signals Systems & Computers, 1-4 Nov. 1998, Vol.2, III, (1998) 959-963. Wang, Q., Ward, R.A New Edge-directed Image Processing, 7-10 Oct. 2001, Vol.3, (2001) 899-902. Bai, Y., Wang, D, On The Comparison of Interpolation Techniques for Robotic Position. IEEE International Conference on Compensation, Systems, Man and Cybernetics, 5-8 Oct., 2003, Vol.4, (2003) 3384-3389. Hadhoud, M. M., Dessouky, M. I., et al, Adaptive Image Interpolation Based on Local Activity Levels. Proceedings of the 20th National Radio Science Conference, 18-20, March, 2003, (2003) C4_1-C4_8. Shenoy, R. G., Parks, T. W., An Optimal Recovery Approach to Interpolation. IEEE Transactions on Signal Processing, Vol.40, (1992) 1987-1996. N. E. Huang et al., The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis. Proc. R. Soc. Lond. A (1998) 454, 903-995. Anna Linderhed, Variable Sampling of the Empirical Mode Decomposition of Two-Dimensional Signals, International Journal of Wavelets, Multiresolution and Information Processing, Vol. 3, No. 3 (September 2005): Special Issue on Frame Theory and Sampling Problems in Time-Frequency Analysis and Wavelet Theory (Part II). Kan Zeng, Ming-Xia He, A Simple Boundary Process Technique for Empirical Mode Decomposition, Geoscience and Remote Sensing Symposium, 2004. IGARSS '04. Proceedings. 2004 IEEE International Vol.6, 2004 Page(s): 4258 - 4261. Tan Tian, Caifang Zhang, An Iterative Hybird Method for Image Interpolation, International Conference on Intelligent Computing (ICIC 2005) (Part I), pp.10-19, 2005.

10040