A combined method for obtaining fringe orientations of ESPI

Optics Communications 273 (2007) 60–66 www.elsevier.com/locate/optcom

A combined method for obtaining fringe orientations of ESPI Xia Yang b

a,*

, Qifeng Yu

a,b

, Sihua Fu

a

a College of Astronautics, National University of Defense Technology, Changsha, Hunan 410073, PR China College of Optical Engineering, National University of Defense Technology, Changsha, Hunan 410073, PR China

Received 31 July 2006; received in revised form 13 December 2006; accepted 15 December 2006

Abstract Fringe orientation can direct the processing of fringe patterns and provide basic information for understanding the fringe patterns. The gradient method is popularly used for its convenience but it is easily affected by noise. The plane-fit method is capable of obtaining precise local fringe orientation despite the high speckle noise but it is strict with the calculating window size. In this paper, a combined method is proposed. This method uses the plane-fit algorithm to get local gradients and uses the gradient method to get the orientation results from these gradients. It can be seen that this new method can get more accurate orientations for ESPI than both the gradient method and the plane-fit method. Ó 2006 Elsevier B.V. All rights reserved. PACS: 42.30.Ms Keywords: ESPI; Fringe orientation; Gradient method; Plane-fit method

1. Introduction Fringe orientation provides important information for understanding or analyzing fringes. It can direct noise filtering [1–4], skeleton extracting [5], or even phase obtaining [6,7]. Yu et al. designed a directional window along the fringe orientation and filtered noise in the directional window. This is called spin filter [1,2]. A similar idea can also be found in InSAR filtering methods [8]. The spin filter uses 8 or 16 discrete orientations and its filtering windows are nearly lines. These limit the efficiency of noise filtering. So Yu et al. developed the spin filter and proposed a contoured window filter [3,9], which designs a contoured window by searching continuous fringe orientations to simulate the fringe contour and then filters noise in the contoured window. They also use fringe orientations to help extracting skeletons [5] and obtaining the whole field phase [6,7]. The basic physical quantity of interest in interferometry is the phase of the observed fringe pattern from which *

Corresponding author. Tel./fax: +86 7314518992. E-mail address: [email protected] (X. Yang).

0030-4018/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2006.12.026

the desired information such as deformations can be derived. In fact, many methods have been proposed to extract the phase from a single fringe image with the help of fringe orientations. Marroquin et al. proposed a quadrature transform method to calculate the phase by estimating local frequencies and orientations [10,11]. To improve this method, they also proposed a fringe-tracking method to get the fringe orientations [12]. Larkin proposed a Hilbert transform method to get the whole filed phase, and one of the key steps is to calculate the fringe orientations [13]. Besides, Fourier transform methods for single interferograms can also benefit from fringe orientations [14]. Various methods of estimating fringe orientations are known from the literature. They include the spin method [1,2], the plane-fit method [3,5,6], the model based method [15,16], the two-dimensional spectral estimation method [17] and the gradient-based method [18–20]. The spin method does well for fringes with low noise. The plane-fit method can get accurate orientations in spite of high noise, but it is strict with the size of the calculating window. The model based method can do well only for some certain

X. Yang et al. / Optics Communications 273 (2007) 60–66

61

types of fringes. The gradient-based method is popularly used for its efficiency, but its accuracy can easily be lowed by the speckle noise of ESPI. In the following parts of this paper it will be shown that our method, which is combined by the plane-fit method and the gradient-based method, is more adaptive than the plane-fit method and can get more accurate results than both methods. 2. Plane-fit method The plane-fit method is proposed here for comparison. More information can be fond in Refs. [3] or [5]. Suppose the fitted plane is gðx; yÞ ¼ a þ bx þ cy

ð1Þ

By least-square fitting, one can get 8 P P P P a 1þb xþc y ¼ fxy > > > x;y x;y x;y x;y > > P 2 P P < P a x þ b x þ c xy ¼ fxy x x;y x;y x;y x;y > > P P P 2 P > > > fxy y : a y þ b xy þ c y ¼ x;y

x;y

x;y

Fig. 1. A real ESPI image.

ð2Þ

x;y

where x and y are coordinates and fxy is the gray value of the point (x, y). Because the selectedParea isPsymmetric Pboth for axis x and axis y, so the items x;y x, x;y y and x;y xy are equal to zero. And in a rectangle window, it is also true that X X x2 ¼ y2 ð3Þ x;y

x;y

So Eq. (2) can be simplified as 8 P P a 1 ¼ fxy > > > x;y x;y > > < P 2 P b x ¼ fxy x x;y x;y > > P 2 P > > > fxy y :c y ¼ x;y

Fig. 2. The orientation result obtained with the plane-fit method.

ð4Þ

x;y

The fringe orientation can be defined as the direction of the normal vector: " P   ogðx; yÞ ogðx; yÞ c x;y fxy y h ¼ arctg ¼ arctg ¼ arctg P oy ox b x;y fxy x x¼y¼0

As shown in Fig. 3, for a local point of O, if the window size is close to the fringe size, the fitted plane is slant (line AB in Fig. 3) and is immune to noise. In this case, good orientation results can be obtained. If the window size is enlarged to double of the fringe size, the fitted plane will be close to horizontal (line CD in Fig. 3) and the results are easily affected by noise. When the window size is enlarged to triple of the fringe size, the fitted plane became slant again (line EF in Fig. 3). This can explain why the plane-fit method is so strict with the size

ð5Þ For a real ESPI image as shown in Fig. 1, the orientation result of the plane-fit method is shown in Fig. 2. Gray 0 represents angle 0 and gray 255 represents angle p. The calculating window size used here is 25  25. From Fig. 2 it can be seen that the plane-fit method can get fine orientations even with high speckle noise. This proves that this method is robust to the high ESPI noise. From Fig. 2 it also can be seen that on ridges and peaks of the fringes, plane-fit method cannot get the right orientations because the fitted planes at these points are horizontal. When the fitted planes are close to horizontal, the orientation results are not good because the normal vector directions of the fitted planes are easily affected by noise.

Fig. 3. The sections of the fitted planes on fringes.

62

X. Yang et al. / Optics Communications 273 (2007) 60–66

of the calculating windows. For local point A or B which is on the ridge and vale of the fringe, the planes fitted there are always close to horizontal as the selected areas are always nearly symmetric. So the plane-fit method cannot get right orientations at ridges and vales of the fringes. 3. Gradient-based method Gradient-based method is popularly used for analyzing fringes [18–20]. Ideally, the fringe direction vector is approximately perpendicular to the intensity gradient vector, so the inner product of these two vectors equals to zero: ðcos h; sin hÞrIðx; yÞ ¼ 0

ð6Þ

Fig. 4. The orientation result obtained with the gradient-based method directly.

where I is the intensity of fringes and h(x, y) is the fringe local tangential direction. So the fringe orientation h can be defined as    oIðx; yÞ oIðx; yÞ hðx; yÞ ¼ arctg ð7Þ oy ox

P ð2I x I y Þ 1 h ¼ arctg P 2 2 ðI y  I 2x Þ

To get more robust results, an average strategy is applied. The notion ‘‘direction’’ is always used to describe angels with a range from 0 to 2p, and the notion ‘‘orientation’’ represents the angle with a range from 0 to p. So opposite directions means the same orientation. What we get from Eq. (7) is the orientation field. The orientations cannot be averaged directly because the same orientations such as 0 and p will annihilate each other. When averaging this field, the local orientation should be doubled so that same orientations such as 0 and p will be the same directions as 0 and 2p, perpendicular orientations such as 0 and p/2 will be opposite directions as 0 and p. In this way, same orientations will enforce each other and perpendicular orientations will annihilate each other. The following equations are used to substitute Eq. (7). The averaging of the doubled orientations can be done in complex domain to avoid the 2p jumps. Suppose a complex data which is

where Ix and Iy are the gradients (derivatives) and can be obtained by operators such as Sobel or Prewitt. The accuracy can be estimated by Cði; jÞ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 1 2 Cði; jÞ ¼ jhðu; vÞ  hði; jÞj N ðu;vÞ2s  d; f ðd ¼ ðh0  h þ 360Þmod360Þ < 180 jhðu; vÞ  hði; jÞj ¼ d  180 otherwise

cðx; yÞ ¼ cos 2hðx; yÞ þ i sin 2hðx; yÞ

ð8Þ

where i is the square root of 1, h(x, y) is the orientations to be averaged and c(x, y) is a complex corresponding to every point in the averaging window. The averaged value of the complex data is P 1X cðx; yÞ 1 X c ¼ cos 2hðx; yÞ þ i sin 2hðx; yÞ ð9Þ ¼ n n n The sums of Eq. (9) are all done in the averaging window. So the orientation can be obtained from the averaged complex value of c: P P 1 sin 2h 1 ð2 sin h cos hÞ h ¼ arctg P ð10Þ ¼ arctg P 2 cos 2h 2 ðcos2 h  sin2 hÞ Substituting (7) into (10), one can get

Ix ¼

oIðx; yÞ ; ox

Iy ¼

oIðx; yÞ oy

ð11Þ ð12Þ

ð13Þ For the same image of Fig. 1, the orientation result of the gradient-based method is shown in Fig. 4. The calculating window size used here is also 25  25. The gradients in Eq. (12) are obtained with Prewitt operate in 3  3 windows. From Fig. 4 it can be seen that with high speckle noise of ESPI, normal gradient-based method cannot get right orientations directly. 4. Combined method As mentioned before in this paper, the plane-fit method is strict with the calculating window size and cannot get right orientations at ridges and vales of the fringes. But the method is robust to high speckle noise. The gradientbased method is very sensitive to noise and cannot get right orientations with the high speckle noise of ESPI directly. So we propose here to combine them to get good orientation results. From the comparison between Eqs. (5) and (7) it can be seen, the plane-fit method computes the gradients too but in a different way. And because the gradients derived by the plane-fit method are more immune to noise, they can be used in Eq. (11) to substitute the derivatives of the gradient-based method. So one can get

X. Yang et al. / Optics Communications 273 (2007) 60–66

63

Fig. 5. The orientation result of Fig. 1 obtained by our method.

P P 1 1 ð2I x I y Þ ð2bcÞ ¼ arctg P 2 h ¼ arctg P 2 2 ðI y  I 2x Þ 2 ðb  c2 Þ where

P x;y fxy y Ix ¼ c ¼ P 2 ; x;y y

P x;y fxy x Iy ¼ b ¼ P 2 x;y x

ð14Þ Fig. 6. The left one is the simulated fringe image with high speckle noise and the right one is its true orientation.

ð15Þ

The b and c in Eq. (15) is the plane-fit parameters of Eq. (1) and can be gotten by Eq. (4). So the steps of our combined method are as follows: (1) For each point in the fringe image, fit a plane by Eq. (1) in the calculating window. (2) Get b and c by Eq. (4) with least-square plane-fit method. (3) For each point, substitute b and c into Eq. (15) and obtain the local orientation by (14) in the averaging window. Fig. 5 is the orientation result of Fig. 1 obtained with our method. The window size is 25  25 for both planefit step and average step. The result of Fig. 5 is better than either the result of the plane-fit method or the result of the gradient-based method. Fig. 5 also shows that the orientations at ridges and vales of the fringes can be correctly estimated.

including the p jumps of orientations. For example, 0 and p are equal orientations, so their differences should be 0 (estimated by Eq. (16)) rather than p. With the plane-fit method, the gradient-based method and the method proposed in this paper, different window sizes are used respectively from 5  5 to about 67  67 with an increment of 2 pixels in every step. In this way not only accuracy can be compared, but also the flexibility to the sizes of the calculating windows. The window size for the plane-fit step of our method is 15  15. Fig. 7 is the result errors of the three methods compared with true orientations. The x-axis is the window size in pixel. The y-axis is the errors calculated by Eq. (16) corresponding to every window size. From Fig. 7 it can be seen that our method has the minimum errors. The plane-fit method is sensitive to window sizes and its result errors are small only when the window sizes are close to integral times of the local fringe width. The gradient method cannot

Gradient-based method Plane-fit method Our method

0.7

5. Experiment results 0.6

The left of Fig. 6 is a simulated fringe image with high speckle noise (40%), and the right of Fig. 6 is the true orientations of the left image. The simulated fringe width is about 28 pixels. To verify the analysis above, three methods are used respectively to obtain orientations from Fig. 6. The orientation errors are calculated by 1X E¼ ð16Þ jsin½hðx; yÞ  hT ðx; yÞj n x;y2s where n is the number of points in area s, hðx; yÞ is the obtained orientation and hT ðx; yÞ is the true orientation shown in the right of Fig. 6. Errors are estimated this way to avoid

0.5 0.4 0.3 0.2 0.1 0 5

9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 Fig. 7. The accuracy comparison of the three methods.

64

X. Yang et al. / Optics Communications 273 (2007) 60–66

get fine results because its gradients are badly affected by high speckle noise. Fig. 8 is the orientation result of our method obtained from the simulated image of Fig. 6. The window size for plane-fit is 15  15, and the window for the complex filter is 25  25. The error computed by Eq. (16) is 0.0182. It should be made clear that the complex filter is equal to the sine–cosine filter essentially, and out method realize the plane-fit method and the complex filter by Eq. (14). Fig. 9 is obtained by the gradient method with a window of 25  25 after filtering by the median filter with a window of 15  15. The error is 0.1607. It can be seen that the median filter combined with the gradient method is not as good as the method proposed in this paper. Fig. 10 is obtained by the

gradient method with a window of 25  25 first and then filtered by 15  15 the sine–cosine filter. It can be seen from Fig. 10 that without good estimations of the gradients, the sine–cosine filter also cannot get good orientation results. And this is the reason why the plane-fit method is chosen in this paper to calculate the gradients, as plane-fit method is robust to noise and can obtain good gradients. Fig. 11 is another low quality ESPI fringe image with high noise. Fig. 12 is the orientation result obtained by our method from Fig. 11. To get more smooth result, the image has been filtered before the calculating of fringe orientations. The window size for plane-fit is 9  9 and the window size for complex filtering is 45  45.

Fig. 8. The orientation image computed by our method.

Fig. 10. The orientation image computed by the gradient method and the sine/cosine filter.

Fig. 9. The orientation image computed by the median filter and the gradient method.

Fig. 11. Another real ESPI fringe image.

X. Yang et al. / Optics Communications 273 (2007) 60–66

65

Fig. 14. The filtering result of Fig. 1 by the contoured window filter.

Fig. 12. The orientation result of Fig. 11.

6. An application of fringe orientation: contoured window filter Based on fringe orientations, the contoured window filter [3,5] can be designed to filter ESPI noise effectively. The contoured window filter filters noise in the contoured windows, and the contoured windows are defined according to the fringe orientations. Suppose the coordinates of the current point is P 0 ðx0 ; y 0 Þ with the local fringe orientation h0. Then one can get two neighboring points P 1 ðx1 ; y 1 Þ and P 1 ðx1 ; y 1 Þ along the orientation and its adverse respectively. Continue on, . . . P n ðxn ; y n Þ and P n ðxn ; y n Þ can be calculated as follows and as shown in Fig. 13a: ( xi ¼ xði1Þ  cos hði1Þ ð17Þ y i ¼ y ði1Þ  sin hði1Þ where 1 6 i 6 n and hi is the local orientation of point ðxi ; y i Þ. Note that when come to the point ðxi ; y i Þ, the next tracing orientation is changed to hi. Tracing points one by one, one can obtain a curve which fits to the equal-phase curve. By widening the curve in both sides the contoured window is determined, as shown in Fig. 13b. Fig. 13c demonstrates that the contoured windows in the right fit well

Fig. 15. The filtering result of Fig. 11 by the contoured window filter.

the equal-phase lines. Because the pixel values are quite homogeneous in the contoured windows, noise can be easily filtered by low pass filers. Fig. 14 is the filtering result of Fig. 1 with the contoured window filter and Fig. 15 is the filtering result of Fig. 11. The filtering results show that with the help of fringe orientations, high speckle noise can be well removed. 7. Conclusions

Fig. 13. The contoured windows obtained by tracing along the fringe orientations.

Fringe orientation is a kind of very important information for fringe processing. But it is not easy to get high accuracy orientations for ESPI when the noise is considerably high. In this paper a new method is proposed to obtain the fringe orientations. This method combines the plane-fit method and the gradient-based method together. The plane-fit method is robust to high speckle noise but it is sensitive to the window size. The gradient-based method is effective and is much adaptive to window size but it fails when the noise is high. So after combination, our method can obtain better orientations than both

66

X. Yang et al. / Optics Communications 273 (2007) 60–66

methods as the experiment results shown. Compared with the gradient-based method, our method is more accurate and more robust to noise. And compared with the planefit method, our method is also more accurate and more adaptive to the window size. The gradients of our proposed method are obtained by the plane-fit method. So the windows for plane-fit should be smaller than the fringe widths to get stable results. Especially when the fringe width varies a lot, the window size must be smaller than the minimum fringe width. Our future work will be put on how to realize this plane-fit step adaptively. References [1] [2] [3] [4]

Qifeng Yu, Appl. Opt. 27 (18) (1988) 3782. Qifeng Yu, X. Liu, K. Andresen, Appl. Opt. 33 (17) (1994) 3705. Qifeng Yu, X. Sun, X. Liu, Appl. Opt. 41 (2002) 2650. G. Kaufmann, A. Davila, D. Kerr, Optical Testing Digest 2 (4) (1997). [5] Qifeng Yu, X. Sun, X. Liu, Opt. Eng. 42 (2003) 68.

[6] Qifeng Yu, Sihua Fu, Xia Yang, Xiangyi Sun, Xiaolin Liu, Opt. Express 12 (2004) 75. . [7] Qifeng Yu, Sihua Fu, Xiaolin Liu, Xia Yang, Xiangyi Sun, Opt. Express 12 (20) (2004) 4980. [8] Jong-Sen Lee, K.P. Papathanassiou, T.L. Ainsworth, M.R. Grunes, A. Reigber, IEEE Trans. Geosci. Remote Sens. 36 (5) (1998). [9] Qifeng Yu, Xia Yang, Sihua Fu, Appl. Opt. 44 (2005) 7050. [10] J.L. Marroquin, M. Servin, R. Rodriguez-Vera, J. Opt. Soc. Am. A 14 (1997) 1742. [11] J.L. Marroquin, R. Rodriguez-Vera, M. Servin, J. Opt. Soc. Am. A 15 (1998) 1536. [12] Juan Antonio Quiroga, Manuel Servin, Francisco Cuevas, J. Opt. Soc. Am. A 19 (8) (2002) 1524. [13] K.G.Larkin,D.J.Bone, M.A. Oldfield,J. Opt. Soc.Am. A 18(2001) 1862. [14] T. Kreis, Akademie-Verlag, Berlin, 1996. [15] P. Vizcaya, L. Gerhardt, Pattern Recognit. 29 (7) (1996) 1221. [16] Jie Zhou, Jinwei Gu, IEEE Trans. Image Process. 13 (6) (2004). [17] C.L. Wilson, G.T. Candela, C.I. Watson, J. Artif. Neural Networks 1 (2) (1994) 203. [18] Xiang Zhou, John P. Baird, John F. Arnold, Appl. Opt. 38 (5) (1999). [19] A.K. Jain, L. Hong, S. Pankanti, R. Bolle, Proc. IEEE 85 (9) (1997) 1365. [20] Asker M. Bazen, Sabih H. Gerez, IEEE Trans. Pattern Anal. Mach. Intell. 24 (7) (2002).