- Email: [email protected]
Comparison of Fourier transform, windowed Fourier transform, and wavelet transform methods for phase calculation at discontinuities in fringe projection profilometry
Optics and Lasers in Engineering 50 (2012) 1152–1160
Contents lists available at SciVerse ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
Comparison of Fourier transform, windowed Fourier transform, and wavelet transform methods for phase calculation at discontinuities in fringe projection profilometry Zonghua Zhang a,n, Zhao Jing a, Zhaohui Wang a, Dengfeng Kuang b a b
School of Mechanical Engineering, Hebei University of Technology, Tianjin 300130, China Institute of Modern Optics, the Key Laboratory of Optical Information Science and Technology of the Education Ministry of China, Nankai University, Tianjin 300071, China
a r t i c l e i n f o
a b s t r a c t
Available online 26 March 2012
Phase demodulation techniques from one fringe pattern have been widely studied because it can measure dynamic objects by capturing single image. These techniques mainly include Fourier transform (FT), windowed Fourier transform (WFT), and wavelet transform (WT). FT has been widely used to demodulate phase information from single deformed fringe pattern on smooth objects. However, for objects having discontinuities and/or large slopes, FT cannot obtain correct phase at the edges because of its global processing. WFT and WT have been applied to nonstationary fringe pattern analysis. Since local fringe information used to extract phase information, WFT and WT are better than FT for phase calculation at discontinuities and/or slopes. In this paper, we discuss the pro and con of the three methods on phase calculation at discontinuities and/or slopes. Simulated and experimental data are tested at edges in order to confirm which method is appropriate to measure objects having discontinuities by using one-frame fringe pattern acquisition method. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Fringe projection profilometry Discontinuities measurement Fourier transform Windowed Fourier transform Wavelet transform
1. Introduction Non-contact optical full-field 3D sensing techniques are active research fields in machine vision, industrial inspection, reverse engineering, physical profiling, biomedical testing, and many other fields [1,2]. Among them, phase calculation-based fringe projection techniques are widely studied because of the advantages of fast acquisition and accurate data retrieval. Two traditional methods for phase calculation are phase-shifting profilometry (PSP) [3] and Fourier transform profilometry (FTP) [4]. PSP method processes the fringe patterns pixel by pixel, but it needs at least three fringe patterns having constant phase shift in between and this technique is sensitive to noise and environment. Therefore, it is hard to measure moving objects by using PSP-based methods. Fourier transform (FT) method has been widely used to demodulate phase information from single deformed fringe pattern on smooth objects, so it has been applied to such fields as fast acquisition 3D shape of unsteady objects. However, for objects having discontinuities and/or large slopes, FT cannot obtain correct phase at the edges because of its global processing.
n
Corresponding author. Tel./fax: þ 86 22 60204189. E-mail address: [email protected] (Z. Zhang).
0143-8166/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.optlaseng.2012.03.004
All frequency components in frequency domain are decided by the entire spatial information, so the higher frequency components at the edges may distribute in the whole wrapped phase map. The obtained phase and then shape data are incorrect because of the frequency leakage. In order to demodulate nonstationary and noisy fringe patterns, recently, windowed Fourier transform (WFT) [5] and wavelet transform (WT) [6] profilometry are proposed. Since they use local fringe information to extract phase information, the two techniques are better than FT for phase calculation at discontinuities and/or slopes. Most recently, Huang et al. compared FT, WFT and WT methods for phase extraction from a single fringe pattern on smooth surface [7]. Although they pointed out that all the three transform-based algorithms perform poorly at edges, there are no further researches investigated regarding the differences among FT, WFT and WT algorithms. Gorthi et al. mainly discussed the effects of non-periodicity of fringe pattern on errors of phase calculation and some qualitative comparisons between FT and WFT at discontinuities [8]. Li et al. investigated different wavelets (1D and 2D) for phase recovery on a hemispherical object [9]. Quan et al. mainly gave theoretical description, simulated and experimental results of phase-shifting, FT, WFT and WT algorithms, and then briefly discussed their performance on phase calculation at discontinuities [10]. However, no one has given the quantitative evaluation and analysis at discontinuities of a single fringe pattern among FT, WFT and WT algorithms.
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
From a fringe pattern, FT, WFT or WT algorithm obtains the wrapped phase map, which needs to be unwrapped for 3D shape conversion. In order to avoid effects of spatial phase unwrapping on shape data at discontinuities, the obtained wrapped phase map by FT, WFT and WT needs to be unwrapped pixel by pixel. Therefore, the optimum three-frequency selection method [11–13] is applied to calculate absolute (or unwrapped) phase map of a fringe pattern. This method uses three fringe patterns pffiffiffiffiffiffi having the optimum fringe numbers of N0, N0 1, and N0 N 0 to calculate the absolute fringe order on a pixel by pixel basis and therefore the technique is compatible with objects possessing discontinuities and/or isolated surfaces. The principle of FT, WFT and WTP is introduced in the following Section. They are used to retrieve the wrapped phase from a single fringe pattern projected on objects having discontinuities. Simulated data having a step are implemented by the three techniques in Section 3 to show which one is suitable to measure discontinuities. Section 4 demonstrates actual experiments on measuring a step by the three techniques. Last Section gives concluding remarks.
2. Principle The three FT, WFT and WT profilometries generally contain the following three steps: (1) Projecting system projects sinusoidal fringe pattern onto the surface of measured objects; (2) From a different viewpoint of the projector, CCD camera acquires the deformed fringe pattern which is modulated by the height of measured object surface to get discrete image information; (3) Calculating the wrapped phase information by FT, WFT or WT algorithm and then the unwrapped phase map by the optimum three-frequency selection method. One sinusoidal fringe pattern projected onto a measured object surface can be represented as f ðx,yÞ ¼ aðx,yÞ þ bðx,yÞcos½jðx,yÞ þ Dfðx,yÞ
2.1. Fourier transform Fourier transform (FT) is a widely used phase demodulation method from one fringe pattern image. In most cases, FT uses fast Fourier transform (FFT) to calculate the wrapped phase data. It includes applying discrete FFT to obtain the spectrum of modulated fringe pattern, extracting the fundamental spectrum, implementing inverse FT, and eventually retrieving the phase information of the measured object. FT and inverse FT can be expressed as Z 1Z 1 Ff ðe, ZÞ ¼ f ðx,yÞexpðjexjZyÞdxdy, ð2Þ
1 4p2
Where Im() and Re() mean imaginary and real part of c(x,y). The obtained phase from Eq. (5) is in the range of p to p, which needs to be unwrapped in order to get the corresponding 3D shape information. Due to the spectral leakage in the neighbourhood of discontinuities and/or the areas of a large surface slope, FT can not give correct phase information [4]. Windowing the fringe pattern before taking FT can reduce the leakage errors [14], therefore, in principle WFT will give more accurate phase information in these areas. Another fault of FT is frequency aliasing, which makes accurate extraction fundamental spectrum difficult. 2.2. Windowed Fourier transform Windowed Fourier transform (WFT) introduces a local window function into FT [5]. The basic idea of WFT is: dividing the whole fringe pattern into a number of local fringe regions with a moving window, obtaining each local fringe pattern’s spectrum by FT; superimposing all the local spectrums to get the whole fringe pattern’s spectrum; extracting the fundamental spectrum; implementing inverse FT, and eventually retrieving the phase information of the fringe pattern. Local spectrums are simpler, so WFT can reduce the effect of frequency aliasing on measurement accuracy. WFT and inverse WFT can be expressed as Z 1Z 1 Sf ðu,v, e, ZÞ ¼ f ðx,yÞg nu,v, e, Z ðx,yÞdxdy, ð6Þ 1
f ðx,yÞ ¼
1 4p2
Z
1 1
Z
1 1
1
Z
1 1
Z
1
Sf ðu,v, e, ZÞ g u,v, e, Z ðx,yÞdedZdudv, 1
ð7Þ g u,v, e, Z ðx,yÞ ¼ gðxu,yvÞexpðjex þjZyÞ,
ð8Þ
gðx,yÞ ¼ exp½x2 =2s2x ðy2 =2sy Þ2 ,
ð9Þ
where ‘n’ denotes the complex conjugate, g(x,y) is a window function. Eq. (9) expresses the Gaussian function. Parameters u and v are the shifting factor of window center in x and y directions, respectively. As u and v increase, the window moves forward to cover the entire domain. sx and sy are the standard deviations of the Gaussian function in x and y directions, respectively, which control the spatial extension of g(x,y). In FT technique, when transforming f(x,y) to Ff(e,Z), the frequency information appears but time information is buried deeply and can hardly be recognized. That is to say, although FT gives frequency values from the spectrum, it is unknown where they occur in the signal. However, from Sf(u,v,e,Z), WFT obtains not only the spectrum components but also where a component appearing in the time domain. So long as the distance between signals is larger than the effective radius of the Gaussian window, the spectrum will not affect each other [5,15]. The obtained phase from WFT is in the range of p to p, which needs to be unwrapped for shape measurement.
1
2.3. Wavelet transform
F I f ðe, ZÞ ¼ BPF½Ff ðe, ZÞ cðx,yÞ ¼
fundamental spectrum of Ff(e,Z). After applying the inverse Fourier transform to the obtained fundamental spectrum by Eq. (4), the wrapped phase can be retrieved as Imcðx,yÞ , ð5Þ fðx,yÞ ¼ arctan Recðx,yÞ
ð1Þ
where f(x,y), a(x,y) and b(x,y) are recorded intensity, background intensity and fringe amplitude, respectively. j(x,y) is the phase of the projected fringe pattern. Df(x,y) is the phase corresponding to the shape of the measured object. In order to demodulate phase Df(x,y) from the captured fringe pattern f(x,y), the principle of FT, WFT and WT is introduced. Because the obtained phase from the three methods is wrapped, the optimum three-frequency selection method is used to independently calculate the absolute phase information at each pixel position.
1
1153
Z
1
1
Z
ð3Þ
1
F I f ðe, ZÞexpðjex þ jZyÞdedZ,
ð4Þ
1
where e and Z are frequency along x and y directions, respectively. BPF[U] is the operator of Band Pass Filter and FIf(e,Z) is the
The principle of the wavelet transform (WT) [6,16] is to perform local filtering operations on a signal with the scaled and shifted versions of a mother wavelet, which is similar to that in FT dividing a signal into its sinusoidal components of various amplitudes and frequencies. Moreover, it is worth mentioning
1154
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
that the wavelet transform has a multi-resolution property in the spatial and frequency domains which reduces the resolution problem inherent in other transforms such as FT and WFT. The chosen mother wavelet in WT should be adapted to the fringe pattern in order to extract correct phase. Examples of analytical wavelets are the complex Morlet, Gaussian and Paul mother wavelets, which are defined by the following equations, respectively. ! 1 x2 ccmor ðxÞ ¼ 2 1=4 expð2pif c xÞexp , ð10Þ 2 2f b ðf b pÞ dðC p expðixÞexpðx2 ÞÞp ccgau ðxÞ ¼ , dxp
ccpau ðxÞ ¼
2n n!ð1ixÞðn þ 1Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , 2p ð2nÞ!=2
ð11Þ
ð12Þ
where fb is the bandwidth parameter and is set here to a value of 1. fc is the mother wavelet centre frequency and it is set here to the value 1.5 or 2. The complex Gaussian wavelet is calculated by taking the pth derivative of the function C p expðixÞexpðx2 Þ. Where Cp is a constant factor which depends on the derivative order p of function C p expðixÞexpðx2 Þ and is computed to normalize the Gaussian wavelet (i.e., :ccgau ðxÞ: ¼ 1). In this paper, p is set to 4. n is the order of the Paul mother wavelet and is here set to a value of 4. Fig. 1 shows 1D complex Morlet, Gaussian and Paul wavelet which will be used in the simulated experiment of Section 3. One-dimensional continuous wavelet transform (1D-CWT) of fringe pattern f(x) can be expressed as Z 1 n W f ða,bÞ ¼ f ðxÞca,b ðxÞdx ¼ /f ðxÞ, ca,b ðxÞS ð13Þ 1
thereamong, 1 xb , ca,b ðxÞ ¼ c a a where ‘n’ denotes the complex conjugate, ca,b(x) is a series of wavelets that obtained through companding and translating the mother wavelet function c(x), a is the scale factor, and b is the shift factor. Eq. (13) describes the inner product of f(x) and ca,b(x)
and the result of Wf(a,b) reflects the similarity of f(x) and ca,b(x). Utilizing analytical wavelets function to do 1D-CWT, the amplitude A(a,b) and phase j(a,b) of the transform can be expressed as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Aða,bÞ ¼ ½imagðW f ða,bÞÞ2 þ ½realðW f ða,bÞÞ2 , ð14Þ imagðW f ða,bÞÞ : realðW f ða,bÞÞ
jða,bÞ ¼ arctan
ð15Þ
The maximum of the WT amplitudes at every line is defined as the ridge of the WT, which can be simply expressed as ridgeðbÞ ¼ maxðAðai ,bÞÞ:
ð16Þ
The phase of the fringe pattern is equal to the phase of its wavelet transform on the ridge of the WT [17]. The obtained phase is wrapped, which needs to be unwrapped for shape calculation as well. 2.4. Optimum three-frequency selection The optimum multi-frequency selection method defines the numbers of projected fringes to be [11]: Nf i ¼ Nf 0 ðNf 0 Þði1Þ=ðn1Þ ,
for,
i ¼ 1,:::,n1
ð17Þ
where Nf0 and Nfi are the maximum number of fringes and the number of fringes in the ith fringe set, respectively, and n is the number of fringe sets used. With three fringe sets, the method is usually referred to as the optimum three-frequency selection method. For example, if Nf0 ¼ 100 and n ¼3, the other two fringe sets have fringe numbers of Nf1 ¼Nf0 1¼ 99 and pffiffiffiffiffiffiffiffi Nf 2 ¼ Nf 0 N f 0 ¼ 90. This method resolves fringe order ambiguity as the beat obtained between Nf0 and Nf1 is a single fringe over the full field of view. The reliability of the obtained fringe order is maximized as fringe order calculation is performed through a geometric series of beat fringes with 1, 10 and 100 fringes. The fringe order calculation for Nf0 ¼ 100 is reliable to 6s (giving a probability of 99.73% of calculating the correct absolute fringe order) providing the phase noise to one standard deviation is better than a certain value of a fringe [12]. The obtained absolute phase data are unwrapped pixel by pixel, so the optimum threefrequency selection method can measure objects with discontinuities and/or isolated surfaces.
Fig. 1. Examples of complex mother wavelets. (a) cmor1-1.5, (b) cmor1-2, (c) cgau4, (d) cpau4.
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
1155
Therefore, in the following simulated and actual experiments, after FT, WFT and WT calculate wrapped phase from a fringe pattern having discontinuities, applying the optimum threefrequency selection method to the wrapped phase data obtains absolute unwrapped phase pixel by pixel. The discontinuities have no effect on the obtained absolute phase data.
3. Comparison by simulated data In order to test FT, WFT and WT methods on measuring objects having discontinuities, simulated data will be used to give quantitative analysis. A step with height of 50 mm was simulated in the software of Matlab. The step is located in the middle of the image with size of 1024 pixel 1024 pixel. Three sinusoidal fringe patterns with the Gaussian noise level of 3.0% were projected onto the step. The direction of fringe patterns is parallel to the step. The three fringe patterns have the optimum fringe numbers of 100, 99 and 90, so the optimum three-frequency selection method can independently calculate the absolute phase pixel by pixel. Fig. 2 illustrates the modulated sinusoidal fringe patterns with fringe numbers of 100 on the step surface. The fringe patterns are vertical and parallel to the step in Fig. 2(a) and (b), respectively. 3.1. 1D-CWT by different wavelets
Fig. 2. The modulated sinusoidal fringe patterns with noise of 3.0% and having fringe numbers of 100 on the step surface. (a) fringe patterns vertical to step and (b) fringe patterns parallel to step.
Due to demodulate the sinusoidal fringe pattern represented by Eq. (1), the wavelets having the similar sinusoidal shape in Fig. 1 were chosen to do 1D-CWT. At first, the wrapped phase data were obtained by doing 1D-CWT to the modulated fringe patterns. Then the optimum three-frequency selection method independently calculates the unwrapped phase at each pixel position. Finally, the unwrapped phase was converted into depth data. In order to check which wavelet is the most suitable to retrieve 3D shape information at discontinuities, the difference between the obtained height and the simulated height was estimated. Fig. 3 shows the average error around the step along middle column direction when fringe patterns are vertical to the step. It can be seen that the error is smaller when using cmor1-1.5 and cmor1-2 wavelet because the 1D wavelet is insensitive to step when the direction of the fringe patterns is vertical to the
Fig. 3. Depth error of simulated step along the middle column direction within the range of pixel position 499 and 525 by using different wavelets for step height 100 mm when fringe patterns vertical to step.
1156
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
step. The main error comes from the noise. Fig. 4 shows the height error along the middle row direction within the range of pixel positions 492 to 532 around the step when fringe patterns are parallel to the step. The average error for the pixel positions 492 to 532 are shown in Table 1 for Morlet1-1.5, Morlet1-2, Paul and Gaussian. It can be shown that at discontinuities Paul wavelet gives the best results, then Gaussian wavelet, and finally Morlet wavelet. The experimental results are the same as the statement in [18] ‘‘the Paul mother wavelet is a wise choice in cases where a fringe pattern with rapid phase variations is being analysed’’. In the following, Paul and Morlet cmor1-1.5 wavelets will be chosen to calculate the wrapped phase data for fringe patterns being parallel and vertical to the step.
PSP, FT, WFT, and WT algorithms. Average height error of the simulated step along middle column direction within the range of pixel positions 492 to 532 by PSP, FT, WFT and WT are 0.07 mm, 5.25 mm, 0.25 mm, 0.56 mm, respectively. It shows that at discontinuities WFT and WT perform better than FT. When fringe patterns are parallel to step, Paul wavelet was used to calculate phase for WT. The height error along the middle row direction is illustrated in Fig. 6 for step height of 50 mm and 80 mm. Table 2 shows the average error for the pixel positions 492 and 532 for the two step heights by PSP, FT, WFT and WT methods. It shows that at discontinuities WFT and WT perform better than FT as well. Because of spectrum leakage, FT gives much wider error data than WFT and WT methods at step position. Although it gives the most accurate depth data at
3.2. Comparison of FT, WFT and WT FT, WFT and WT together with the optimum three-frequency selection method were utilized to obtain the height of the simulated step. The obtained height of the step was compared to that of the simulated step. When fringe patterns are vertical to step, Morlet (cmor1-1.5) wavelet was used to calculate phase for WT. Fig. 5 illustrates the reconstructed depth error of the simulated step with height 100 mm along the middle column direction within the range of pixel position 492 to 532 by using
Table 1 Average height error of the simulated step along middle row direction within the range of pixel positions 492 and 532 for step height of 50 mm and 80 mm by using different mother wavelets (Unit mm). Step height (mm)
Morlet1-1.5
Morlet1-2
Paul
Gaussian
50 80
5.00 13.22
5.37 16.31
2.96 4.35
4.03 7.03
Fig. 4. Depth error along the middle row direction within the range of pixel position 492 and 532 of the simulated step by using different mother wavelets for two different height steps when fringe patterns parallel to step. (a) 50 mm and (b) 80 mm.
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
1157
Fig. 5. Depth error of simulated step along the middle column direction within the range of pixel position 492 and 532 by using PSP, FT, WFT, and WT algorithms for step height 100 mm when fringe patterns vertical to step.
Fig. 6. Depth error along the middle row direction within the range of pixel position 492 and 532 by using PSP, FT, WFT, and WT for two different height steps when fringe patterns parallel to step. (a) 50 mm and (b) 80 mm.
steps, PSP needs more than three fringe patterns to obtain the wrapped phase data. From the simulated experimental results, it is concluded that the error of the reconstructed data at
discontinuities is small when the fringe patterns are vertical to steps. Therefore, in actual measurements, it is best to project fringe patterns vertical to steps. In the following experiments, the
1158
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
projected fringe patterns are vertical to steps to calculate the phase data. It should be pointed out that the integration limits and the window extension must be noticed when using WFT. They may affect the quality of the wrapped phase seriously [19]. After numerous experiments, we get the following facts: (1) The period of the fringe pattern needs to be estimated and then the center of the integration limits e0 is calculated roughly by [20] 2p e0 ¼ , Tx
ð18Þ
where Tx is the average period calculated from the fringe pattern. A step means more high frequencies, thererfore upper integration limit should have a larger value in order to correctly calculate phase information at discontinuities. In
Table 2 Average height error of the simulated step along middle row direction within the range of pixel positions 492 and 532 for step height of 50 mm and 80 mm by PSP, FT, WFT and WT methods. Step height (mm)
PSP
FT
WFT
WT
50 80
0.08 0.08
3.86 6.27
1.34 2.04
2.96 4.35
the simulated step data, when there are 100, 99 and 90 fringes on two flat surfaces, in principle the center integration limit e0 are 0.6136, 0.6075, 0.5522, respectively. However, the upper limit of 2.5 was chosen in order to include more frequencies at the step position. (2) By comparing the different values of window extension s, it needs to be a little larger than the period of carrier fringe pattern. 4. Experimental comparison 4.1. Hardware system The hardware setup comprises a portable DLP (Digital Light Processing) video projector, a 3-CCD color camera with IEEE 1394 port and a personal computer (PC), as illustrated in Fig. 7. The projector is from BenQ (Model CP270) with one-chip digital micro-mirror device (DMD) and a resolution of up to 1024 pixel 768 pixel (XGA). The 3-CCD camera from Hitachi (Model HVF22F) has a resolution of 1360 pixel 1024 pixel. In order to completely avoid the crosstalk between color channels, we use a single colour channel. Three green images with different spatial frequency is generated in the PC and projected onto an object surface, for example a manufactured step, by using an uneven fringe projection method [21]. The camera captures the three fringe pattern images and saves them into the computer for post processing.
4.2. Results
Fig. 7. The hardware setup of the 3D imaging system including a DLP projector, a color 3CCD camera and a personal computer.
A manufactured step having known distance was measured as the target object. Three fringe pattern images with the optimum fringe numbers of 100, 99 and 90 were successively generated in the computer and projected by the DLP projector onto the step’s surface, as demonstrated in Fig. 7. The 3CCD camera captured the three deformed fringe pattern images from another view of point, as illustrated in Fig. 8. Applying FT, WFT, and WT together with the optimum threefrequency selection method to the captured fringe pattern images obtain the unwrapped phase map. Because PSP can accurately calculate phase in the absence of noise, the phase data by fourstep PSP were used as a standard baseline. The unwrapped phase data were converted to depth data by using the developed calibration method [22]. The obtained depth data by FT, WFT and WT methods were compared to that of PSP. Fig. 9 illustrates the reconstructed depth data along one column direction by using PS, FT, WFT and WT algorithms. In order to clearly show the difference among the four algorithms, 5 mm up and down shift
Fig. 8. Three fringe pattern images of a manufactured step having different fringe numbers. (a) 90, (b) 99, and (c) 100.
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
1159
Fig. 9. Reconstructed depth data along one column direction by using PS, FT, WFT and WT methods. (a) profile of the obtained depth data (5 mm up and down shift in order to clearly show the difference), (b) the partial enlargement of (a).
Fig. 10. The whole 3D reconstructed depth map by using four different algorithms. (a) PSP, (b) FT, (c) WFT and (d) WT.
Fig. 11. Difference of the reconstructed depth map between PSP and FT, WFT, WT. (a) FT, (b) WFT and (c) WT.
was added to the profile in Fig. 9(a). Fig. 10 illustrates the whole 3D reconstructed depth map by using the four algorithms and Fig. 11 shows the difference of the reconstructed depth data between PSP and FT, WFT, WT. The average error between PSP and FT, WFT, WT for the pixel positions of 460 to 471 around one step is 2.91 mm, 0.89 mm and 1.17 mm, respectively. The experimental results show that at discontinuities of the step, WFT gives reliable phase data information, and then WT
algorithm, while FT can not obtain correct phase information. The algorithm of FT, WFT and WT were written in the software of Matlab and implemented on a personal computer. The configuration of the computer is Intel(R) Core(TM) 2 Duo P7350 CPU with 2.00 GHZ and 2 GB memory. The processing time of FT, WFT, and WT is 13 s, 348 s and 120 s, respectively.
1160
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
5. Conclusions
References
In this paper, three one-frame fringe pattern demodulation methods of FT, WFT, and WT are discussed to calculate phase information at discontinuities. For measuring objects having discontinuities, WFT is the most appropriate transform method to retrieve wrapped phase from one fringe pattern. The reason is WFT uses a Gaussian filter window at each pixel position to calculate phase information. However, the computation cost is much more expensive than FT method. Because of spectrum leakage at discontinuities, FT causes large error and can not give correct phase information. However, the processing time is much less than WFT and WT methods. Therefore, for measuring objects having smooth surface by one-frame fringe pattern projection technique, FT is the first choice; while for objects having discontinuities, the projected fringe patterns should be vertical to the step and then WFT can give the best results with the cost of longer computation time than FT. With the development of new software algorithms and advent of hardware such as graphics processing unit (GPU), all the oneframe transform-based methods can greatly reduce processing time in the near future. Therefore, they are promising methods in real time 3D imaging system and have many applications in the fields of multi-media, medical imaging, on-line quality inspection, human-computer interaction (HCI) and security industries. Further study on more accurately extracting absolute phase information and then converting to 3D shape data are under the way.
[1] Chen F, Brown GM, Song M. Overview of three-dimensional shape measurement using optical methods. Opt Eng 2000;39:10–22. [2] Blais F. Review of 20 years of range sensor development. J Electron Imaging 2004;13:231–40. [3] Phillion DW. General methods for generating phase-shifting interferometry algorithm. Appl Opt 1997;36:8098–115. [4] Su XY, Chen WJ. Fourier transform profilometry: a review. Opt Lasers Eng 2001;35:263–84. [5] Qian KM. Two-dimensional windowed Fourier transform for fringe pattern analysis: principles, applications and implementations. Opt Lasers Eng 2007;45:304–17. [6] A.Z.A. Abid, Fringe pattern analysis using Wavelet Transforms, Ph. D. Thesis, Liverpool John Moores University (2008). [7] Huang L, Qian KM, Pan B, Asundi AK. Comparison of Fourier transform, windowed Fourier transform, and wavelet transform methods for phase extraction from a single fringe pattern in fringe projection profilometry. Opt Lasers Eng 2010;48:141–8. [8] Gorthi SS, Lolla KR. Wavelet transform analysis of truncated fringe patterns in 3D surface profilometry. Proc SPIE 2005;5856:265. [9] Li SK, Su XY, Chen WJ. Spatial carrier fringe pattern phase demodulation by use of a two-dimensional real wavelet. Appl Opt 2009;48:6893–906. [10] Quan C, Chen W, Tay CJ. Phase-retrieval techniques in fringe-projection profilometry. Opt Lasers Eng 2010;48:235–43. [11] Zhang ZH, Towers CE, Towers DP. Time efficient colour fringe projection system for 3D shape and colour using optimum 3-frequency interferometry. Opt Express 2006;14:6444–55. [12] Towers CE, Towers DP, Jones JDC. Generalized frequency selection in multifrequency interferometry. Opt Lett 2004;29:1348–50. [13] Towers CE, Towers DP, Jones JDC. Absolute fringe order calculation using optimised multi-frequency selection in full-field porfilometry. Opt Lasers Eng 2005;43:788–800. [14] Berryman F, Pynsent P, Cubillo J. The effect of windowing in Fourier transform profilometry applied to noisy images. Opt Lasers Eng 2004;41:815–25. [15] Qian KM, Soon SH. Two-dimensional windowed Fourier frames for noise reduction in fringe pattern analysis. Opt Eng 2005;44:075601. [16] Zhong JG, Weng JW. Spatial carrier-fringe pattern analysis by means of wavelet transform: wavelet transform profilometry. Appl Opt 2004;43:4993–8. [17] Zhong JG, Weng JW. Phase retrieval of optical fringe patterns from the ridge of a wavelet transform. Opt Lett 2005;30:2560–2. [18] Gdeisat MA, Abid A, Burton DR, Lalor MJ, Lilley F, Moore C, Qudeisat M. Spatial and temporal carrier fringe pattern demodulation using the onedimensional continuous wavelet transform: Recent process, challenges, and suggested developments. Opt Lasers Eng 2009;47:1348–61. [19] Qian KM. On window size selection in the windowed Fourier ridges algorithm. Opt Lasers Eng 2008;46:1186–92. [20] Cheng P, Hu JS, Zhang GF. Deformation measurements of dragonfly’s wings in free flight by using windowed Fourier transform. Opt Lasers Eng 2008;46:157–61. [21] Zhang ZH, Towers CE, Towers DP. Uneven fringe projection for efficient calibration in high resolution 3D shape metrology. Appl Opt 2007;46:6113–9. [22] Zhang ZH, Ma HY, Zhang SX, Guo T, Towers CE, Towers DP. Simple calibration of a phase-based 3D imaging system based on uneven fringe projection. Opt Lett 2011;36:627–9.
Acknowledgments The authors would like to thank the National Natural Science Foundation of China (61171048, 10904076), Program for New Century Excellent Talents in University (NCET), the Key Project of Chinese Ministry of Education (No: 211016), Specialized Research Fund for the Doctoral Program of Higher Education (’’SRFDP’’) (No: 20111317120002), the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry (NO: 20101561), Research Project supported by Hebei Education Department (No: ZD2010121), and Scientific Research Foundation for the Returned Overseas Scholars, Hebei Province.
Contents lists available at SciVerse ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
Comparison of Fourier transform, windowed Fourier transform, and wavelet transform methods for phase calculation at discontinuities in fringe projection profilometry Zonghua Zhang a,n, Zhao Jing a, Zhaohui Wang a, Dengfeng Kuang b a b
School of Mechanical Engineering, Hebei University of Technology, Tianjin 300130, China Institute of Modern Optics, the Key Laboratory of Optical Information Science and Technology of the Education Ministry of China, Nankai University, Tianjin 300071, China
a r t i c l e i n f o
a b s t r a c t
Available online 26 March 2012
Phase demodulation techniques from one fringe pattern have been widely studied because it can measure dynamic objects by capturing single image. These techniques mainly include Fourier transform (FT), windowed Fourier transform (WFT), and wavelet transform (WT). FT has been widely used to demodulate phase information from single deformed fringe pattern on smooth objects. However, for objects having discontinuities and/or large slopes, FT cannot obtain correct phase at the edges because of its global processing. WFT and WT have been applied to nonstationary fringe pattern analysis. Since local fringe information used to extract phase information, WFT and WT are better than FT for phase calculation at discontinuities and/or slopes. In this paper, we discuss the pro and con of the three methods on phase calculation at discontinuities and/or slopes. Simulated and experimental data are tested at edges in order to confirm which method is appropriate to measure objects having discontinuities by using one-frame fringe pattern acquisition method. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Fringe projection profilometry Discontinuities measurement Fourier transform Windowed Fourier transform Wavelet transform
1. Introduction Non-contact optical full-field 3D sensing techniques are active research fields in machine vision, industrial inspection, reverse engineering, physical profiling, biomedical testing, and many other fields [1,2]. Among them, phase calculation-based fringe projection techniques are widely studied because of the advantages of fast acquisition and accurate data retrieval. Two traditional methods for phase calculation are phase-shifting profilometry (PSP) [3] and Fourier transform profilometry (FTP) [4]. PSP method processes the fringe patterns pixel by pixel, but it needs at least three fringe patterns having constant phase shift in between and this technique is sensitive to noise and environment. Therefore, it is hard to measure moving objects by using PSP-based methods. Fourier transform (FT) method has been widely used to demodulate phase information from single deformed fringe pattern on smooth objects, so it has been applied to such fields as fast acquisition 3D shape of unsteady objects. However, for objects having discontinuities and/or large slopes, FT cannot obtain correct phase at the edges because of its global processing.
n
Corresponding author. Tel./fax: þ 86 22 60204189. E-mail address: [email protected] (Z. Zhang).
0143-8166/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.optlaseng.2012.03.004
All frequency components in frequency domain are decided by the entire spatial information, so the higher frequency components at the edges may distribute in the whole wrapped phase map. The obtained phase and then shape data are incorrect because of the frequency leakage. In order to demodulate nonstationary and noisy fringe patterns, recently, windowed Fourier transform (WFT) [5] and wavelet transform (WT) [6] profilometry are proposed. Since they use local fringe information to extract phase information, the two techniques are better than FT for phase calculation at discontinuities and/or slopes. Most recently, Huang et al. compared FT, WFT and WT methods for phase extraction from a single fringe pattern on smooth surface [7]. Although they pointed out that all the three transform-based algorithms perform poorly at edges, there are no further researches investigated regarding the differences among FT, WFT and WT algorithms. Gorthi et al. mainly discussed the effects of non-periodicity of fringe pattern on errors of phase calculation and some qualitative comparisons between FT and WFT at discontinuities [8]. Li et al. investigated different wavelets (1D and 2D) for phase recovery on a hemispherical object [9]. Quan et al. mainly gave theoretical description, simulated and experimental results of phase-shifting, FT, WFT and WT algorithms, and then briefly discussed their performance on phase calculation at discontinuities [10]. However, no one has given the quantitative evaluation and analysis at discontinuities of a single fringe pattern among FT, WFT and WT algorithms.
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
From a fringe pattern, FT, WFT or WT algorithm obtains the wrapped phase map, which needs to be unwrapped for 3D shape conversion. In order to avoid effects of spatial phase unwrapping on shape data at discontinuities, the obtained wrapped phase map by FT, WFT and WT needs to be unwrapped pixel by pixel. Therefore, the optimum three-frequency selection method [11–13] is applied to calculate absolute (or unwrapped) phase map of a fringe pattern. This method uses three fringe patterns pffiffiffiffiffiffi having the optimum fringe numbers of N0, N0 1, and N0 N 0 to calculate the absolute fringe order on a pixel by pixel basis and therefore the technique is compatible with objects possessing discontinuities and/or isolated surfaces. The principle of FT, WFT and WTP is introduced in the following Section. They are used to retrieve the wrapped phase from a single fringe pattern projected on objects having discontinuities. Simulated data having a step are implemented by the three techniques in Section 3 to show which one is suitable to measure discontinuities. Section 4 demonstrates actual experiments on measuring a step by the three techniques. Last Section gives concluding remarks.
2. Principle The three FT, WFT and WT profilometries generally contain the following three steps: (1) Projecting system projects sinusoidal fringe pattern onto the surface of measured objects; (2) From a different viewpoint of the projector, CCD camera acquires the deformed fringe pattern which is modulated by the height of measured object surface to get discrete image information; (3) Calculating the wrapped phase information by FT, WFT or WT algorithm and then the unwrapped phase map by the optimum three-frequency selection method. One sinusoidal fringe pattern projected onto a measured object surface can be represented as f ðx,yÞ ¼ aðx,yÞ þ bðx,yÞcos½jðx,yÞ þ Dfðx,yÞ
2.1. Fourier transform Fourier transform (FT) is a widely used phase demodulation method from one fringe pattern image. In most cases, FT uses fast Fourier transform (FFT) to calculate the wrapped phase data. It includes applying discrete FFT to obtain the spectrum of modulated fringe pattern, extracting the fundamental spectrum, implementing inverse FT, and eventually retrieving the phase information of the measured object. FT and inverse FT can be expressed as Z 1Z 1 Ff ðe, ZÞ ¼ f ðx,yÞexpðjexjZyÞdxdy, ð2Þ
1 4p2
Where Im() and Re() mean imaginary and real part of c(x,y). The obtained phase from Eq. (5) is in the range of p to p, which needs to be unwrapped in order to get the corresponding 3D shape information. Due to the spectral leakage in the neighbourhood of discontinuities and/or the areas of a large surface slope, FT can not give correct phase information [4]. Windowing the fringe pattern before taking FT can reduce the leakage errors [14], therefore, in principle WFT will give more accurate phase information in these areas. Another fault of FT is frequency aliasing, which makes accurate extraction fundamental spectrum difficult. 2.2. Windowed Fourier transform Windowed Fourier transform (WFT) introduces a local window function into FT [5]. The basic idea of WFT is: dividing the whole fringe pattern into a number of local fringe regions with a moving window, obtaining each local fringe pattern’s spectrum by FT; superimposing all the local spectrums to get the whole fringe pattern’s spectrum; extracting the fundamental spectrum; implementing inverse FT, and eventually retrieving the phase information of the fringe pattern. Local spectrums are simpler, so WFT can reduce the effect of frequency aliasing on measurement accuracy. WFT and inverse WFT can be expressed as Z 1Z 1 Sf ðu,v, e, ZÞ ¼ f ðx,yÞg nu,v, e, Z ðx,yÞdxdy, ð6Þ 1
f ðx,yÞ ¼
1 4p2
Z
1 1
Z
1 1
1
Z
1 1
Z
1
Sf ðu,v, e, ZÞ g u,v, e, Z ðx,yÞdedZdudv, 1
ð7Þ g u,v, e, Z ðx,yÞ ¼ gðxu,yvÞexpðjex þjZyÞ,
ð8Þ
gðx,yÞ ¼ exp½x2 =2s2x ðy2 =2sy Þ2 ,
ð9Þ
where ‘n’ denotes the complex conjugate, g(x,y) is a window function. Eq. (9) expresses the Gaussian function. Parameters u and v are the shifting factor of window center in x and y directions, respectively. As u and v increase, the window moves forward to cover the entire domain. sx and sy are the standard deviations of the Gaussian function in x and y directions, respectively, which control the spatial extension of g(x,y). In FT technique, when transforming f(x,y) to Ff(e,Z), the frequency information appears but time information is buried deeply and can hardly be recognized. That is to say, although FT gives frequency values from the spectrum, it is unknown where they occur in the signal. However, from Sf(u,v,e,Z), WFT obtains not only the spectrum components but also where a component appearing in the time domain. So long as the distance between signals is larger than the effective radius of the Gaussian window, the spectrum will not affect each other [5,15]. The obtained phase from WFT is in the range of p to p, which needs to be unwrapped for shape measurement.
1
2.3. Wavelet transform
F I f ðe, ZÞ ¼ BPF½Ff ðe, ZÞ cðx,yÞ ¼
fundamental spectrum of Ff(e,Z). After applying the inverse Fourier transform to the obtained fundamental spectrum by Eq. (4), the wrapped phase can be retrieved as Imcðx,yÞ , ð5Þ fðx,yÞ ¼ arctan Recðx,yÞ
ð1Þ
where f(x,y), a(x,y) and b(x,y) are recorded intensity, background intensity and fringe amplitude, respectively. j(x,y) is the phase of the projected fringe pattern. Df(x,y) is the phase corresponding to the shape of the measured object. In order to demodulate phase Df(x,y) from the captured fringe pattern f(x,y), the principle of FT, WFT and WT is introduced. Because the obtained phase from the three methods is wrapped, the optimum three-frequency selection method is used to independently calculate the absolute phase information at each pixel position.
1
1153
Z
1
1
Z
ð3Þ
1
F I f ðe, ZÞexpðjex þ jZyÞdedZ,
ð4Þ
1
where e and Z are frequency along x and y directions, respectively. BPF[U] is the operator of Band Pass Filter and FIf(e,Z) is the
The principle of the wavelet transform (WT) [6,16] is to perform local filtering operations on a signal with the scaled and shifted versions of a mother wavelet, which is similar to that in FT dividing a signal into its sinusoidal components of various amplitudes and frequencies. Moreover, it is worth mentioning
1154
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
that the wavelet transform has a multi-resolution property in the spatial and frequency domains which reduces the resolution problem inherent in other transforms such as FT and WFT. The chosen mother wavelet in WT should be adapted to the fringe pattern in order to extract correct phase. Examples of analytical wavelets are the complex Morlet, Gaussian and Paul mother wavelets, which are defined by the following equations, respectively. ! 1 x2 ccmor ðxÞ ¼ 2 1=4 expð2pif c xÞexp , ð10Þ 2 2f b ðf b pÞ dðC p expðixÞexpðx2 ÞÞp ccgau ðxÞ ¼ , dxp
ccpau ðxÞ ¼
2n n!ð1ixÞðn þ 1Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , 2p ð2nÞ!=2
ð11Þ
ð12Þ
where fb is the bandwidth parameter and is set here to a value of 1. fc is the mother wavelet centre frequency and it is set here to the value 1.5 or 2. The complex Gaussian wavelet is calculated by taking the pth derivative of the function C p expðixÞexpðx2 Þ. Where Cp is a constant factor which depends on the derivative order p of function C p expðixÞexpðx2 Þ and is computed to normalize the Gaussian wavelet (i.e., :ccgau ðxÞ: ¼ 1). In this paper, p is set to 4. n is the order of the Paul mother wavelet and is here set to a value of 4. Fig. 1 shows 1D complex Morlet, Gaussian and Paul wavelet which will be used in the simulated experiment of Section 3. One-dimensional continuous wavelet transform (1D-CWT) of fringe pattern f(x) can be expressed as Z 1 n W f ða,bÞ ¼ f ðxÞca,b ðxÞdx ¼ /f ðxÞ, ca,b ðxÞS ð13Þ 1
thereamong, 1 xb , ca,b ðxÞ ¼ c a a where ‘n’ denotes the complex conjugate, ca,b(x) is a series of wavelets that obtained through companding and translating the mother wavelet function c(x), a is the scale factor, and b is the shift factor. Eq. (13) describes the inner product of f(x) and ca,b(x)
and the result of Wf(a,b) reflects the similarity of f(x) and ca,b(x). Utilizing analytical wavelets function to do 1D-CWT, the amplitude A(a,b) and phase j(a,b) of the transform can be expressed as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Aða,bÞ ¼ ½imagðW f ða,bÞÞ2 þ ½realðW f ða,bÞÞ2 , ð14Þ imagðW f ða,bÞÞ : realðW f ða,bÞÞ
jða,bÞ ¼ arctan
ð15Þ
The maximum of the WT amplitudes at every line is defined as the ridge of the WT, which can be simply expressed as ridgeðbÞ ¼ maxðAðai ,bÞÞ:
ð16Þ
The phase of the fringe pattern is equal to the phase of its wavelet transform on the ridge of the WT [17]. The obtained phase is wrapped, which needs to be unwrapped for shape calculation as well. 2.4. Optimum three-frequency selection The optimum multi-frequency selection method defines the numbers of projected fringes to be [11]: Nf i ¼ Nf 0 ðNf 0 Þði1Þ=ðn1Þ ,
for,
i ¼ 1,:::,n1
ð17Þ
where Nf0 and Nfi are the maximum number of fringes and the number of fringes in the ith fringe set, respectively, and n is the number of fringe sets used. With three fringe sets, the method is usually referred to as the optimum three-frequency selection method. For example, if Nf0 ¼ 100 and n ¼3, the other two fringe sets have fringe numbers of Nf1 ¼Nf0 1¼ 99 and pffiffiffiffiffiffiffiffi Nf 2 ¼ Nf 0 N f 0 ¼ 90. This method resolves fringe order ambiguity as the beat obtained between Nf0 and Nf1 is a single fringe over the full field of view. The reliability of the obtained fringe order is maximized as fringe order calculation is performed through a geometric series of beat fringes with 1, 10 and 100 fringes. The fringe order calculation for Nf0 ¼ 100 is reliable to 6s (giving a probability of 99.73% of calculating the correct absolute fringe order) providing the phase noise to one standard deviation is better than a certain value of a fringe [12]. The obtained absolute phase data are unwrapped pixel by pixel, so the optimum threefrequency selection method can measure objects with discontinuities and/or isolated surfaces.
Fig. 1. Examples of complex mother wavelets. (a) cmor1-1.5, (b) cmor1-2, (c) cgau4, (d) cpau4.
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
1155
Therefore, in the following simulated and actual experiments, after FT, WFT and WT calculate wrapped phase from a fringe pattern having discontinuities, applying the optimum threefrequency selection method to the wrapped phase data obtains absolute unwrapped phase pixel by pixel. The discontinuities have no effect on the obtained absolute phase data.
3. Comparison by simulated data In order to test FT, WFT and WT methods on measuring objects having discontinuities, simulated data will be used to give quantitative analysis. A step with height of 50 mm was simulated in the software of Matlab. The step is located in the middle of the image with size of 1024 pixel 1024 pixel. Three sinusoidal fringe patterns with the Gaussian noise level of 3.0% were projected onto the step. The direction of fringe patterns is parallel to the step. The three fringe patterns have the optimum fringe numbers of 100, 99 and 90, so the optimum three-frequency selection method can independently calculate the absolute phase pixel by pixel. Fig. 2 illustrates the modulated sinusoidal fringe patterns with fringe numbers of 100 on the step surface. The fringe patterns are vertical and parallel to the step in Fig. 2(a) and (b), respectively. 3.1. 1D-CWT by different wavelets
Fig. 2. The modulated sinusoidal fringe patterns with noise of 3.0% and having fringe numbers of 100 on the step surface. (a) fringe patterns vertical to step and (b) fringe patterns parallel to step.
Due to demodulate the sinusoidal fringe pattern represented by Eq. (1), the wavelets having the similar sinusoidal shape in Fig. 1 were chosen to do 1D-CWT. At first, the wrapped phase data were obtained by doing 1D-CWT to the modulated fringe patterns. Then the optimum three-frequency selection method independently calculates the unwrapped phase at each pixel position. Finally, the unwrapped phase was converted into depth data. In order to check which wavelet is the most suitable to retrieve 3D shape information at discontinuities, the difference between the obtained height and the simulated height was estimated. Fig. 3 shows the average error around the step along middle column direction when fringe patterns are vertical to the step. It can be seen that the error is smaller when using cmor1-1.5 and cmor1-2 wavelet because the 1D wavelet is insensitive to step when the direction of the fringe patterns is vertical to the
Fig. 3. Depth error of simulated step along the middle column direction within the range of pixel position 499 and 525 by using different wavelets for step height 100 mm when fringe patterns vertical to step.
1156
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
step. The main error comes from the noise. Fig. 4 shows the height error along the middle row direction within the range of pixel positions 492 to 532 around the step when fringe patterns are parallel to the step. The average error for the pixel positions 492 to 532 are shown in Table 1 for Morlet1-1.5, Morlet1-2, Paul and Gaussian. It can be shown that at discontinuities Paul wavelet gives the best results, then Gaussian wavelet, and finally Morlet wavelet. The experimental results are the same as the statement in [18] ‘‘the Paul mother wavelet is a wise choice in cases where a fringe pattern with rapid phase variations is being analysed’’. In the following, Paul and Morlet cmor1-1.5 wavelets will be chosen to calculate the wrapped phase data for fringe patterns being parallel and vertical to the step.
PSP, FT, WFT, and WT algorithms. Average height error of the simulated step along middle column direction within the range of pixel positions 492 to 532 by PSP, FT, WFT and WT are 0.07 mm, 5.25 mm, 0.25 mm, 0.56 mm, respectively. It shows that at discontinuities WFT and WT perform better than FT. When fringe patterns are parallel to step, Paul wavelet was used to calculate phase for WT. The height error along the middle row direction is illustrated in Fig. 6 for step height of 50 mm and 80 mm. Table 2 shows the average error for the pixel positions 492 and 532 for the two step heights by PSP, FT, WFT and WT methods. It shows that at discontinuities WFT and WT perform better than FT as well. Because of spectrum leakage, FT gives much wider error data than WFT and WT methods at step position. Although it gives the most accurate depth data at
3.2. Comparison of FT, WFT and WT FT, WFT and WT together with the optimum three-frequency selection method were utilized to obtain the height of the simulated step. The obtained height of the step was compared to that of the simulated step. When fringe patterns are vertical to step, Morlet (cmor1-1.5) wavelet was used to calculate phase for WT. Fig. 5 illustrates the reconstructed depth error of the simulated step with height 100 mm along the middle column direction within the range of pixel position 492 to 532 by using
Table 1 Average height error of the simulated step along middle row direction within the range of pixel positions 492 and 532 for step height of 50 mm and 80 mm by using different mother wavelets (Unit mm). Step height (mm)
Morlet1-1.5
Morlet1-2
Paul
Gaussian
50 80
5.00 13.22
5.37 16.31
2.96 4.35
4.03 7.03
Fig. 4. Depth error along the middle row direction within the range of pixel position 492 and 532 of the simulated step by using different mother wavelets for two different height steps when fringe patterns parallel to step. (a) 50 mm and (b) 80 mm.
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
1157
Fig. 5. Depth error of simulated step along the middle column direction within the range of pixel position 492 and 532 by using PSP, FT, WFT, and WT algorithms for step height 100 mm when fringe patterns vertical to step.
Fig. 6. Depth error along the middle row direction within the range of pixel position 492 and 532 by using PSP, FT, WFT, and WT for two different height steps when fringe patterns parallel to step. (a) 50 mm and (b) 80 mm.
steps, PSP needs more than three fringe patterns to obtain the wrapped phase data. From the simulated experimental results, it is concluded that the error of the reconstructed data at
discontinuities is small when the fringe patterns are vertical to steps. Therefore, in actual measurements, it is best to project fringe patterns vertical to steps. In the following experiments, the
1158
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
projected fringe patterns are vertical to steps to calculate the phase data. It should be pointed out that the integration limits and the window extension must be noticed when using WFT. They may affect the quality of the wrapped phase seriously [19]. After numerous experiments, we get the following facts: (1) The period of the fringe pattern needs to be estimated and then the center of the integration limits e0 is calculated roughly by [20] 2p e0 ¼ , Tx
ð18Þ
where Tx is the average period calculated from the fringe pattern. A step means more high frequencies, thererfore upper integration limit should have a larger value in order to correctly calculate phase information at discontinuities. In
Table 2 Average height error of the simulated step along middle row direction within the range of pixel positions 492 and 532 for step height of 50 mm and 80 mm by PSP, FT, WFT and WT methods. Step height (mm)
PSP
FT
WFT
WT
50 80
0.08 0.08
3.86 6.27
1.34 2.04
2.96 4.35
the simulated step data, when there are 100, 99 and 90 fringes on two flat surfaces, in principle the center integration limit e0 are 0.6136, 0.6075, 0.5522, respectively. However, the upper limit of 2.5 was chosen in order to include more frequencies at the step position. (2) By comparing the different values of window extension s, it needs to be a little larger than the period of carrier fringe pattern. 4. Experimental comparison 4.1. Hardware system The hardware setup comprises a portable DLP (Digital Light Processing) video projector, a 3-CCD color camera with IEEE 1394 port and a personal computer (PC), as illustrated in Fig. 7. The projector is from BenQ (Model CP270) with one-chip digital micro-mirror device (DMD) and a resolution of up to 1024 pixel 768 pixel (XGA). The 3-CCD camera from Hitachi (Model HVF22F) has a resolution of 1360 pixel 1024 pixel. In order to completely avoid the crosstalk between color channels, we use a single colour channel. Three green images with different spatial frequency is generated in the PC and projected onto an object surface, for example a manufactured step, by using an uneven fringe projection method [21]. The camera captures the three fringe pattern images and saves them into the computer for post processing.
4.2. Results
Fig. 7. The hardware setup of the 3D imaging system including a DLP projector, a color 3CCD camera and a personal computer.
A manufactured step having known distance was measured as the target object. Three fringe pattern images with the optimum fringe numbers of 100, 99 and 90 were successively generated in the computer and projected by the DLP projector onto the step’s surface, as demonstrated in Fig. 7. The 3CCD camera captured the three deformed fringe pattern images from another view of point, as illustrated in Fig. 8. Applying FT, WFT, and WT together with the optimum threefrequency selection method to the captured fringe pattern images obtain the unwrapped phase map. Because PSP can accurately calculate phase in the absence of noise, the phase data by fourstep PSP were used as a standard baseline. The unwrapped phase data were converted to depth data by using the developed calibration method [22]. The obtained depth data by FT, WFT and WT methods were compared to that of PSP. Fig. 9 illustrates the reconstructed depth data along one column direction by using PS, FT, WFT and WT algorithms. In order to clearly show the difference among the four algorithms, 5 mm up and down shift
Fig. 8. Three fringe pattern images of a manufactured step having different fringe numbers. (a) 90, (b) 99, and (c) 100.
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
1159
Fig. 9. Reconstructed depth data along one column direction by using PS, FT, WFT and WT methods. (a) profile of the obtained depth data (5 mm up and down shift in order to clearly show the difference), (b) the partial enlargement of (a).
Fig. 10. The whole 3D reconstructed depth map by using four different algorithms. (a) PSP, (b) FT, (c) WFT and (d) WT.
Fig. 11. Difference of the reconstructed depth map between PSP and FT, WFT, WT. (a) FT, (b) WFT and (c) WT.
was added to the profile in Fig. 9(a). Fig. 10 illustrates the whole 3D reconstructed depth map by using the four algorithms and Fig. 11 shows the difference of the reconstructed depth data between PSP and FT, WFT, WT. The average error between PSP and FT, WFT, WT for the pixel positions of 460 to 471 around one step is 2.91 mm, 0.89 mm and 1.17 mm, respectively. The experimental results show that at discontinuities of the step, WFT gives reliable phase data information, and then WT
algorithm, while FT can not obtain correct phase information. The algorithm of FT, WFT and WT were written in the software of Matlab and implemented on a personal computer. The configuration of the computer is Intel(R) Core(TM) 2 Duo P7350 CPU with 2.00 GHZ and 2 GB memory. The processing time of FT, WFT, and WT is 13 s, 348 s and 120 s, respectively.
1160
Z. Zhang et al. / Optics and Lasers in Engineering 50 (2012) 1152–1160
5. Conclusions
References
In this paper, three one-frame fringe pattern demodulation methods of FT, WFT, and WT are discussed to calculate phase information at discontinuities. For measuring objects having discontinuities, WFT is the most appropriate transform method to retrieve wrapped phase from one fringe pattern. The reason is WFT uses a Gaussian filter window at each pixel position to calculate phase information. However, the computation cost is much more expensive than FT method. Because of spectrum leakage at discontinuities, FT causes large error and can not give correct phase information. However, the processing time is much less than WFT and WT methods. Therefore, for measuring objects having smooth surface by one-frame fringe pattern projection technique, FT is the first choice; while for objects having discontinuities, the projected fringe patterns should be vertical to the step and then WFT can give the best results with the cost of longer computation time than FT. With the development of new software algorithms and advent of hardware such as graphics processing unit (GPU), all the oneframe transform-based methods can greatly reduce processing time in the near future. Therefore, they are promising methods in real time 3D imaging system and have many applications in the fields of multi-media, medical imaging, on-line quality inspection, human-computer interaction (HCI) and security industries. Further study on more accurately extracting absolute phase information and then converting to 3D shape data are under the way.
[1] Chen F, Brown GM, Song M. Overview of three-dimensional shape measurement using optical methods. Opt Eng 2000;39:10–22. [2] Blais F. Review of 20 years of range sensor development. J Electron Imaging 2004;13:231–40. [3] Phillion DW. General methods for generating phase-shifting interferometry algorithm. Appl Opt 1997;36:8098–115. [4] Su XY, Chen WJ. Fourier transform profilometry: a review. Opt Lasers Eng 2001;35:263–84. [5] Qian KM. Two-dimensional windowed Fourier transform for fringe pattern analysis: principles, applications and implementations. Opt Lasers Eng 2007;45:304–17. [6] A.Z.A. Abid, Fringe pattern analysis using Wavelet Transforms, Ph. D. Thesis, Liverpool John Moores University (2008). [7] Huang L, Qian KM, Pan B, Asundi AK. Comparison of Fourier transform, windowed Fourier transform, and wavelet transform methods for phase extraction from a single fringe pattern in fringe projection profilometry. Opt Lasers Eng 2010;48:141–8. [8] Gorthi SS, Lolla KR. Wavelet transform analysis of truncated fringe patterns in 3D surface profilometry. Proc SPIE 2005;5856:265. [9] Li SK, Su XY, Chen WJ. Spatial carrier fringe pattern phase demodulation by use of a two-dimensional real wavelet. Appl Opt 2009;48:6893–906. [10] Quan C, Chen W, Tay CJ. Phase-retrieval techniques in fringe-projection profilometry. Opt Lasers Eng 2010;48:235–43. [11] Zhang ZH, Towers CE, Towers DP. Time efficient colour fringe projection system for 3D shape and colour using optimum 3-frequency interferometry. Opt Express 2006;14:6444–55. [12] Towers CE, Towers DP, Jones JDC. Generalized frequency selection in multifrequency interferometry. Opt Lett 2004;29:1348–50. [13] Towers CE, Towers DP, Jones JDC. Absolute fringe order calculation using optimised multi-frequency selection in full-field porfilometry. Opt Lasers Eng 2005;43:788–800. [14] Berryman F, Pynsent P, Cubillo J. The effect of windowing in Fourier transform profilometry applied to noisy images. Opt Lasers Eng 2004;41:815–25. [15] Qian KM, Soon SH. Two-dimensional windowed Fourier frames for noise reduction in fringe pattern analysis. Opt Eng 2005;44:075601. [16] Zhong JG, Weng JW. Spatial carrier-fringe pattern analysis by means of wavelet transform: wavelet transform profilometry. Appl Opt 2004;43:4993–8. [17] Zhong JG, Weng JW. Phase retrieval of optical fringe patterns from the ridge of a wavelet transform. Opt Lett 2005;30:2560–2. [18] Gdeisat MA, Abid A, Burton DR, Lalor MJ, Lilley F, Moore C, Qudeisat M. Spatial and temporal carrier fringe pattern demodulation using the onedimensional continuous wavelet transform: Recent process, challenges, and suggested developments. Opt Lasers Eng 2009;47:1348–61. [19] Qian KM. On window size selection in the windowed Fourier ridges algorithm. Opt Lasers Eng 2008;46:1186–92. [20] Cheng P, Hu JS, Zhang GF. Deformation measurements of dragonfly’s wings in free flight by using windowed Fourier transform. Opt Lasers Eng 2008;46:157–61. [21] Zhang ZH, Towers CE, Towers DP. Uneven fringe projection for efficient calibration in high resolution 3D shape metrology. Appl Opt 2007;46:6113–9. [22] Zhang ZH, Ma HY, Zhang SX, Guo T, Towers CE, Towers DP. Simple calibration of a phase-based 3D imaging system based on uneven fringe projection. Opt Lett 2011;36:627–9.
Acknowledgments The authors would like to thank the National Natural Science Foundation of China (61171048, 10904076), Program for New Century Excellent Talents in University (NCET), the Key Project of Chinese Ministry of Education (No: 211016), Specialized Research Fund for the Doctoral Program of Higher Education (’’SRFDP’’) (No: 20111317120002), the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry (NO: 20101561), Research Project supported by Hebei Education Department (No: ZD2010121), and Scientific Research Foundation for the Returned Overseas Scholars, Hebei Province.