Multi-frequency fringe projection profilometry based on wavelet transform

Chao Jiang; Shuhai Jia; Jun Dong; Qin Lian; Dichen Li

doi:10.1364/OE.24.011323

1. Introduction

Because of its great potential in scientific research and industrial applications, much research has been devoted to three-dimensional (3D)-profiling of objects [1]. With non-contact measurements performed at a large distance, fringe projection profilometry has become one of the alternative techniques [2,3]. Typically, the main technical difficulty is calculating the phase maps from the objects and reference planes. Distinct from phase-measuring profilometry, which needs four or even more images at a time to achieve enough phase information, Fourier transform profilometry (FTP) is more appropriate for dynamic measurements because it only needs 1 or 2 images to calculate phase information [4]. However, as FTP is based on filtering requiring only a single spectrum corresponding to the fundamental frequency component, spectral components may overlap making complete extraction of the first harmonic impossible [5]. The technique of using phase-locked loops (PLLs) in the fringe analysis has developed a lot owing to its well performance [6]. However, this performance depends on the choice of the closed loop gain which is experience-based and lies on the bandwidth of the fringe pattern. In recent years, applying wavelet transform (WT) theory in fringe projection profilometry has been the focus of a great deal of research [7–10]. Because of dilation and translation characteristics of the daughter wavelets, WTs bring time-frequency multi-resolution analytics so that the signal can be well represented in both time-space and frequency domains. Compared with FTP, WT is more suited in actual measurements to phase information retrieval from a non-stationary signal.

Another technical difficulty arising with fringe projection profilometry is the unwrapping algorithm of the phase maps. Because the majority of the unwrapping methods are path-dependent, errors may accumulate and these methods no longer apply in profiling objects with large separated discontinuities as the phase of these discontinuities (multiple of 2 $π$ ) exceeds the range of the unwrapping method (− $π$ to $π$ ) [11–13]. Although the introduction of PLLs replace the conventional phase unwrapping process, it requires that the phase change introduced by the object is smooth and continuous [6]. Meanwhile, a novel multi-frequency fringe projection profilometry was proposed [14,15] in which a peak searching algorithm is adopted to profile an object merely based on the phase information of each point from phase maps for multiple frequency. By using this algorithm, the conventional unwrapping process is avoided and largely discontinuous objects can be measured. However, in regard to the retrieval of the phase map, the adopted FTP cannot achieve a satisfactory performance in actual experiments compared with results obtained using WT.

In this paper, we present a WT-based fringe projection profilometry. A 1D wavelet transform is chosen to obtain the phase maps for each frequency and a peak searching algorithm is applied to obtain the height of the object. This is achieved by independently constructing and solving an object function based on the phases at different frequency for a single point on the object. Owing to the usage of the wavelet transform instead of FTP, the extraction of the phase information is precise enough. Furthermore, by utilizing the peak searching algorithm, the object height of a single measuring point is calculated directly and independently rather than on the basis of the absolute phase from the wrapped phase of the neighbored measuring points. Hence, our method is capable of rapidly profiling the object with large discontinuities.

2. Principle of the multi-frequency color-marked fringe projection profilometry

In a typical fringe projection profilometry system (Fig. 1), a fringe pattern with multiple frequency is directly projected onto the object to be measured using a DMD projector; a CCD camera records the reflected fringe pattern. A fixed plane supports the object and provides a reference plane when calculating the true height of the object.

Fig. 1 System layout of multi-frequency fringe projection profilometry on WT.

Download Full Size | PDF

2.1 Obtaining the wrapped phase map using the WT

The optical geometry associated with fringe projection profilometry (Fig. 2) establishes a coordinate system with the X and Y axes aligned along the vertical and horizontal directions, respectively, of a reference plane. Given $x$ and $y$ as the coordinates of a sampled point on the reference plane, then $h (x, y)$ denotes the height of the object’s surface above the plane. We set $d$ as the distance between the optic centers of the CCD camera and the DMD projector, and $l_{0}$ the perpendicular distance of the optic center of the CCD camera to the plane.

Fig. 2 The optical geometry of the fringe projection profilometry.

Download Full Size | PDF

In the process of measurement, the original projected fringe pattern and the modulated reflected fringe pattern are written respectively as:

I (x, y) = a (x, y) + b (x, y) \cos [2 π f_{0} x + φ_{I} (x, y)]

O (x, y) = a (x, y) + b (x, y) \cos [2 π f_{0} x + φ_{O} (x, y)] .

Here

a (x, y)

is the intensity of the background and

b (x, y)

is the modulation intensity;

f_{0}

represents the fundamental frequency of the fringe pattern, and

φ_{O} (x, y)

and

φ_{I} (x, y)

represent, respectively, the phase modulated by the object and the phase without the object on the reference plane. For simplicity, the projected fringes are aligned parallel to the

y

axis and therefore there is no modulation in this direction.

To obtain the phase map, we chose the complex Morlet wavelet defined as

ψ (x) = π^{1 / 4} \exp (i c x) \exp (- x^{2} / 2) .

Here

c

is the fixed spatial frequency of the Morlet wavelet, the value of which is 5 or 6 to achieve better transform performance. This mother wavelet gives a better localization in both spatial and frequency domains as well as a better optimal window shape [9, 10]. The relative daughter wavelets

ψ_{b, x}

are defined by

ψ_{b, x} (x) = \frac{1}{s} ψ (\frac{x - b}{s}),

where

s

and

b

are the dilation and translation factors. By changing their values, the center frequency of the Morlet wavelet decreases gradually while moving the wavelet along the time axis. In this manner, the signal is thus being processed self-adaptively. Because the fringe pattern is applied to the 1D continuous wavelet transform (CWT) algorithm column by column, the result of the WT for a column of the fringe pattern

I (x)

can be expressed as

W (s, b) = \int_{- \infty}^{\infty} ψ^{*}_{b, x} (x) I (x) d x .

This result is an element of a 2D complex array, which provides the time and frequency information of the fringe pattern calculated.

The argument of $W (s, b)$ , $φ (x, b)$ , is calculated using

φ (s, b) = \tan^{- 1} [\frac{i m a g (W (s, b))}{r e a l (W (s, b))}],

where

i m a g (W (s, b))

and

r e a l (W (s, b))

are the imaginary and real parts of

W (s, b)

. By exploiting the dilation and translation of the Morlet wavelet, the position of the fringe pattern

I (x)

in

φ (x, b)

can be calculated knowing the position of the ridge in

W (s, b)

. This ridge is formed by the loci of the maximum for each frequency, which occur when the frequency of the wavelet is close to the frequency of the calculated fringe pattern.

Figure 3 illustrates the process of extracting the phase map of the fringe pattern using the complex Morlet 1D CWT. The measured object is a simulated semi-sphere (Fig. 3(a)) that produces the reflected fringe pattern (Fig. 3(b)) with a frequency of 0.03 fringes per pixel; column 200 of this fringe pattern is indicated by a black line along which the intensity is plotted (Fig. 3(c)). After applying the 1D CWT, the time and frequency information for column 200 is obtained (Fig. 3(d)). Indicating large values in the image, the red and yellow regions combine to form the ridge. The frequency $f$ (vertical axis) of the ridge matches well with that of the transformed fringe pattern. Figure 3(e) illustrates the phase result obtained using Eq. (6). As the maximum value of the modulus for each column in Fig. 3(d) can be directly determined, its corresponding phase can be obtained from Fig. 3(e), shown as a black line. The resulting phase map of the fringe pattern with all columns is shown in Fig. 3(f).

Fig. 3 (a) The semi sphere simulated in MATLAB. (b) The simulated reflected fringe pattern. (c) The relationship between the intensity and pixel position in column 200. (d) Time and frequency information of the sinusoidal fringe in column 200. (e) The phase result of the sinusoidal fringe in column 200 by using Eq. (6). (f) The phase map of the fringe pattern using complex Morlet 1D CWT.

Download Full Size | PDF

Hence, with the above method, the wrapped phase $φ_{I} (x, y)$ and $φ_{O} (x, y)$ can be obtained and the phase difference as modulated by the object,

d φ (x, y) = φ_{O} (x, y) - φ_{I} (x, y),

can be calculated. From the geometrical setup (Fig. 2),

d φ (x, y)

can be expressed and simplified as [15]

d φ (x, y) = 2 π f_{0} d \times H (x, y) / l_{0} .

Consequently, the wrapped phase maps of the reference plane and the object measured are obtained by applying the WT.

2.2 Peak searching algorithm for calculating object height

Considering the fact that the modulus can be taken into account if the periodicity of a sinusoid is used, a function, which gives a peak signal for the height of the object, is formed by simply adding the sinusoids together and can be written as [14]

S (h) = \frac{1}{K} \sum_{k = 1}^{K} \cos [C_{k} h - d φ (x, y)],

where

K

is the number of the frequency used in composing

S (h)

and

h

is a parameter that is varied over a certain range in finding the real height of the object. The coefficient

C_{k}

is written

C_{k} = 2 π f_{k} d / l_{0},

in which

f_{k}

is the spatial frequency of the fringe pattern. We therefore write

S (h) = \frac{1}{K} \sum_{k = 1}^{K} \cos [\frac{2 π f_{k} d}{l_{0}} h - d φ (x, y)] .

Since only the unit-valued cosines and measured phase information are used to find the real height of the object, the position of the peak from Eq. (11) is phase dependent. With the help of the complex notation as well as the magnitude, only the envelop remains, which enable us to accurately locate the maximum of $S (h)$ . Therefore, we rewrite and simplify the function as:

S (h) = \frac{1}{K} | \sum_{k = 1}^{K} \exp {[\frac{2 π f_{k} d}{l_{0}} (h - H)] \times i} | .

Here

H

is the real height of the object. Note that

S (h)

is merely related to the phase modulated by the object height at the measuring point

(x, y)

.

S (h)

attains its maximum and object height is directly and independently calculated when the parameter

h

equals to the real height of the object

H

without calculating the absolute phase according to the wrapped phase of the neighbored points. Hence, the position of the highest peak represents the real height of the object and the calculating process of the object height is independent and avoids the need for unwrapping.

If the difference between successive spatial frequency of the fringe patterns are equal, i.e.

f_{k} = f_{k - 1} + m,

for all 1<k≤K, Eq. (12) simplifies [15] giving

S (h) = \frac{1}{K} \times \frac{\sin (K α)}{\sin (α)},

where

α

is assumed as

α = \frac{2 π d}{l_{0}} \times \frac{m}{2} (h - H) .

Equation (14) features maxima for

α = \pm n π, n = 1, 2, 3 \dots

. This means that

S (h)

has maxima when

h = H \pm (n l_{0} / d m)

. The curve of Eq. (14), plotted in Fig. 4, is known as the peak-searching curve.

Fig. 4 Illustrative curve of the peak signal S(h) of Eq. (14).

Download Full Size | PDF

The distance between the highest peaks (free height) is calculated using [15]

F r e e H e i g h t = \frac{l_{0}}{d m} .

From Fig. 4, we note that within the range of free height, there is only one main maximum, thus guaranteeing the uniqueness of the result while searching the height of the object without encountering other highest peaks. The free height in the peak searching algorithm determines the range of our measuring system under a certain depth of field.

Because of the geometry (Fig. 2), the resolution of our measuring system is written as

r e s o l u t i o n = \frac{ϕ \times l_{0}}{2 π f_{0} \times d},

where

ϕ

is the minimum phase value that can be measured. Under the resolution limit of the CCD camera, the greater the value of

f_{0}

, the finer the resolution.

3. Experimental results

In the experiment, fringe patterns with multiple frequency that were generated by LabVIEW were projected onto the object using a digital micromirror device (DMD) projector while the reflected fringe patterns were recorded using a CCD camera (CM-030GE, JAI Ltd., Shanghai, China). Because of the property of the WT, only two images of the object and the reference plane are needed for each frequency. To ensure the best measurement result, 11 fringe patterns were generated with frequency that were chosen within the range from 1.22 to 1.62 fringes/cm spaced 0.04 fringes/cm apart.

The object to be profiled in the experiment was a work piece with a step and a sloping surface (Fig. 5). Using vernier calipers, the height of the higher surface was measured at 7.21 mm whereas that of the lower surface was 3.52 mm.

Fig. 5 Work piece profiled in the experiment.

Download Full Size | PDF

Figure 6(a) presents a typical image of the reflected fringe pattern for a certain frequency. Using a complex Morlet 1D CWT, the wrapped phase map obtained of the object at the same frequency is shown in Fig. 6(b).

Fig. 6 (a) Image of the reflected fringe pattern with a certain frequency. (b) Wrapped phase map of the object calculated using the complex Morlet 1D CWT.

Download Full Size | PDF

The results for the work piece measured in the experiment are illustrated in Fig. 7. Figure 7(a) gives a gray image of the object; Fig. 7(b) is the reconstructed 3D shape. The average height of the higher step of the reconstructed object is 7.2517 cm whereas the real height is 7.21 cm. The average measured height of the lower step is 3.5390 cm whereas the real height of the lower step is 3.52 cm. In Fig. 8(b), a greater error clearly exists on the higher step of the reconstructed object; one of the possible explanations is interference from the exposure and shadowing effects. The error can be reduced if another angle of incidence of light is adopted or the intensity of light is decreased.

Fig. 7 (a) Gray-scale image of the object measured. (b) Reconstructed 3D shape of the work piece measured.

Download Full Size | PDF

Fig. 8 (a) Reconstructed 3D shape of the work piece calculated using the 5-step phase shifting method. (b) Reconstructed 3D shape of the work piece calculated using the FTP method. (c) Reconstructed 3D shape of the work piece calculated using the Mexican hat 1D CWT. (d) Reconstructed 3D shape of the work piece calculated using the complex Morlet 1D CWT.

Download Full Size | PDF

To highlight the advantages of multi-frequency fringe projection profilometry based on WT, the measurement results obtained with these methods were compared with that calculated using FTP [15]. Distinct from the FTP and the WT methods, no filtering process is needed but more than one image is required in applying the phase shifting method. Hence, the wrapped phase, calculated by the 5-step phase shifting method, is sufficiently accurate and the corresponding reconstruction result can be viewed as providing a standard measurement value to evaluate the measurement results based on the WT and the FTP method. Besides, owing to its well performance in the reconstruction of the object detail and discontinuous area, the real wavelet, 1-D Mexican hat wavelet, is also widely used in the fringe projection profilometry [16]. Hence, the object was also reconstructed by using the 1-D Mexican hat wavelet.

Figures 8(a)–8(d) illustrate the reconstructed 3D shape of the measured work piece using the 5-step phase shifting method, the FTP method, the Mexican hat wavelet 1-D CWT method and the complex Morlet 1D CWT method, respectively.

Figure 9 provides corresponding cross-sectional drawings of the work piece using these four methods where Fig. 9(a) is the results at low noise area while Fig. 9(b) is the results at high noise area. Note that, when the noise is considerable (exposure area), the reconstructed result using complex Morlet wavelet is closer to that using phase shifting method when compared with 1-D Mexican hat wavelet. However, when the collected fringe patterns are of good quality (non-exposure area), 1-D Mexican hat wavelet has a better reconstruction performance than the complex Morlet wavelet. Both of the two wavelets perform better than the FTP method.

Fig. 9 (a) Comparison profiles of a work piece calculated using four distinct methods at low noise area. (b) Comparison profiles of a work piece calculated using four distinct methods at high noise area.

Download Full Size | PDF

Table 1 is the error statistics of the results calculated using FTP method, 1-D Mexican hat wavelet (Mexh) and 1-D complex Morlet wavelet (Cmor) with respect to the standard measurement result using phase shifting method. Maximum deviation (MD) and standard deviation (SD) were calculated, respectively. It can be seen that the result calculated by Mexican hat wavelet has a maximum fluctuation while complex Morlet wavelet has the minimum one. Hence it can be concluded that the multi-frequency fringe projection profilometry based on WT eliminates noise and offers higher measurement accuracy for objects with large-depth discontinuities while comparing with the one based on FTP method and Mexican hat wavelet.

Table 1. The error statistics of the reconstructed results using FTP, Mexh and Cmor

View Table

4. Conclusion

A multi-frequency fringe projection profilometry with WT was described. The wrapped phase maps of the reference plane and the object at multiple frequency are calculated using the complex Morlet 1D CWT. The peak searching algorithm reconstructs the shape of the object measured from the calculated phase information at different frequency. The final measurement result obtained using the WT was compared with that from the conventional FTP method and Mexican hat wavelet, and exhibited lower noise and higher measurement accuracy with similar measuring times, confirming the reliability of the method. Moreover, by using the WT and the peak searching algorithm, our technique has great potential in providing fast 3D profiling of objects with large depth discontinuities.

Acknowledgments

This project is supported by National Natural Science Foundation of China (Grant No. 51575437, and U1510114), the National High Technology Research and Development Program of China (Grant No. 2015AA020303) and Key Science and Technology Program of Shaanxi Province of China (Grant No. 2014K07-02).

References and links

1. F. Chen, G. Brown, and M. Song, “Overview of three-dimensional shape measurement using optical methods,” Opt. Eng. 39(1), 10–222 (2000). [CrossRef]

2. J. Zhong and Y. Zhang, “Absolute phase-measurement technique based on number theory in multifrequency grating projection profilometry,” Appl. Opt. 40(4), 492–500 (2001). [CrossRef] [PubMed]

3. S. S. Gorthi and P. Rastogi, “Fringe projection techniques: whither we are,” Opt. Lasers Eng. 48(2), 133–140 (2010). [CrossRef]

4. X. Y. Su and W. Chen, “Fourier transform profilometry: a review,” Opt. Lasers Eng. 35(5), 263–284 (2001). [CrossRef]

5. J. Li, X. Y. Su, and L. Y. Guo, “Improved Fourier transform profilometry for the automatic measurement of three-dimensional object shapes,” Opt. Eng. 29(12), 1439–1444 (1990). [CrossRef]

6. D. Ganotra, J. Joseph, and K. Singh, “Second- and first-order phase-locked loops in fringe profilometry and application of neural networks for phase-to-depth conversion,” Opt. Commun. 217(1-6), 85–96 (2003). [CrossRef]

7. J. Zhong and J. Weng, “Spatial carrier-fringe pattern analysis by means of wavelet transform: wavelet transform profilometry,” Appl. Opt. 43(26), 4993–4998 (2004). [CrossRef] [PubMed]

8. J. Zhong and J. Weng, “Phase retrieval of optical fringe patterns from the ridge of a wavelet transform,” Opt. Lett. 30(19), 2560–2562 (2005). [CrossRef] [PubMed]

9. P. Tomassini, A. Giulietti, L. A. Gizzi, M. Galimberti, D. Giulietti, M. Borghesi, and O. Willi, “Analyzing laser plasma interferograms with a continuous wavelet transform ridge extraction technique: the method,” Appl. Opt. 40(35), 6561–6568 (2001). [CrossRef] [PubMed]

10. A. Durson, S. Ozder, and N. Ecevit, “Continuous wavelet transform analysis of projected fringe patterns,” Meas. Sci. Technol. 15(9), 1768–1772 (2004). [CrossRef]

11. L. Xiong and S. Jia, “Phase-error analysis and elimination for nonsinusoidal waveforms in Hilbert transform digital-fringe projection profilometry,” Opt. Lett. 34(15), 2363–2365 (2009). [CrossRef] [PubMed]

12. C. A. Hobson, H. T. Atkinson, and F. Lilley, “The application of digital filtering to phase recovery when surface contouring using fringe projection techniques,” Opt. Lasers Eng. 27(4), 355–368 (1997). [CrossRef]

13. X. Y. Su and W. Chen, “Reliability-guided phase unwrapping algorithm: a review,” Opt. Lasers Eng. 42(3), 245–261 (2004). [CrossRef]

14. Y. Xu, S. H. Jia, X. Luo, J. Yang, and Y. Zhang, “Multi-frequency projected fringe profilometry for measuring objects with large depth discontinuities,” Opt. Commun. 288, 27–30 (2013). [CrossRef]

15. Y. Xu, S. Jia, Q. Bao, H. Chen, and J. Yang, “Recovery of absolute height from wrapped phase maps for fringe projection profilometry,” Opt. Express 22(14), 16819–16828 (2014). [CrossRef] [PubMed]

16. X. Zhou and H. Zhao, “Three-dimensional profilometry based on mexican hat wavelet transform,” Acta Opt. Sin. 29(1), 197–202 (2009). [CrossRef]

	TP	Mexh	Cmor
MD/cm	7.8500	10	3.8750
SD/cm	0.3547	0.7359	0.3111

Multi-frequency fringe projection profilometry based on wavelet transform

Abstract

1. Introduction

2. Principle of the multi-frequency color-marked fringe projection profilometry

2.1 Obtaining the wrapped phase map using the WT

2.2 Peak searching algorithm for calculating object height

3. Experimental results

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (9)

Tables (1)

Equations (17)

Optics Express