Recovery of absolute height from wrapped phase maps for fringe projection profilometry

Yong Xu; Shuhai Jia; Qingchen Bao; Hualing Chen; Jia Yang

doi:10.1364/OE.22.016819

1. Introduction

Three-dimensional (3D) shape measurement techniques have been developed to measure the profile of objects for many years [1]. In both scientific research and industrial applications, such as solid modeling and automated manufacturing, non-contact and full-field optical surface profilometry is becoming increasingly important. Different methods have been investigated and developed, among which those based on fringe projection and triangulation have attracted a lot of attention and have been applied widely in practical applications because of their superior performance and simple implementation [2]. In a simple fringe projection system, a digital micro-mirror device (DMD) projector is usually used to project structured fringe patterns onto an object and a charge coupled device (CCD) camera is used to record the distorted fringes. The relative phase map of the object can be calculated with different methods such as phase-measuring profilometry (PMP) [3–5], Fourier transform profilometry (FTP) [6]. Compared with PMP and other method used to acquire the relative phase map of an object, FTP needs only one or two images of the deformed fringe pattern, which makes it possible for real-time data processing [7]. During the calculation, the involved arctangent operation makes the value of the relative phase map range from – π to π. Thus, the relative phase map is wrapped and spatial phase unwrapping algorithms must be used to make the relative phase map continuous. However, if the object itself has large-depth discontinuities (corresponding phase of multiple of 2π), it will be difficult to unwrap the phase correctly [8]. Many factors can cause errors in phase unwrapping, such as local shadows or low fringe modulation, irregular surface brightness, fringe discontinuities and under sampling [9]. Moreover, many of the phase unwrapping algorithms are path-dependent, and in this case, the phase errors will propagate through the image. The nonlinear intensity response of the DMD projector and the CCD camera will also cause phase errors and therefore measurement errors [10]. Intensive studies have been devoted to developing reliable phase unwrapping algorithms, and a review on this topic can be found in [9].

In 1993, J. M. Huntley and H. Saldner proposed a temporal phase unwrapping algorithm for measuring objects with surface discontinuities [11–13]. In this method, multiple images are recorded at fringe patterns with different frequencies (from low to high) in time line, and the phase of each pixel is unwrapped along the time axis. The phase map calculated by the fringe patterns with the lowest frequency ranges from –π to π, while the phase map with lower fringe frequencies is used to guide the phase unwrapping process at a higher fringe frequency. In this way, the phase of each pixel is calculated by the same pixel in different phase maps with different fringe frequencies, and phase errors do not propagate through the whole image. Although the temporal phase unwrapping method can be used in measuring objects with large depth discontinuities, it needs strict control of fringe parameters. To carry out temporal phase unwrapping, the first group of fringe patterns must contain only one fringe in the full-field to make the first phase map unwrapped, and the real phase difference of the same point between two adjacent groups of fringe patterns must be less than 2π. Moreover, fringe projection profilometry techniques which use multi-frequency fringe patterns have been researched actively [14] and [15]. In these multi-frequency profilometry techniques, two or more fringe patterns with different frequencies are used in order to get the wrapped phase map of an object. In [14] and [15], the need of phase unwrapping procedure still remains [16]. reports a profilometry using four-step phase-shift method, but these four images are used to calculate the local illumination image function, not wrapped phase maps. It is not a phase-measuring profilometry, but an improved line scanning method. Heterodyne principle is also used in multi-frequency fringe projection profilometry such as [14] and [17]. In order to obtain an absolute phase map, multi-frequency fringe projection profilometry with heterodyne principle needs to control the frequencies used in the measurement. The frequency have to be chosen carefully in order to avoid unwrapping failures caused by phase errors [17]. Overall, it is still a challenge to measure complex shaped objects with large depth discontinuities by fringe projection profilometry.

In this paper, we propose a novel multi-frequency fringe projection profilometry with Fourier transform. Fourier transform is used to obtain the wrapped phase map of the object and a peak searching algorithm is utilized to calculate the height of object. In this novel method, the height of each point on the object is measured independently and the phase unwrapping procedure is avoided, enabling the measurement of objects with large depth discontinuities.

2. Principle and method of multi-frequency fringe projection profilometry

In a typical fringe projection profilometry system (Fig. 1), we use a computer to generate fringe patterns with multiple frequencies and a DMD projector to project the patterns onto the object. A CCD camera is used to collect the deformed fringe patterns and send the images to the computer.

Fig. 1 System layout of fringe projection profilometry.

Download Full Size | PDF

2.1 Obtain the wrapped phase map by Fourier transform

The cross section of our system geometry is shown in Fig. 2. X axis and Y axis form the reference plane. $h (x, y)$ is the height of the object. After we place the object on the reference plane, the deformed fringe pattern observed through a CCD camera can be written as:

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

When

h (x, y) = 0

, the deformed image can be expressed by:

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

Here

O (x, y)

represents the fringe pattern modulated by the height of the object and

I (x, y)

represents the original fringe pattern projected onto the reference plane.

a (x, y)

is the background intensity and

b (x, y)

is the modulation intensity.

f_{0}

is the fundamental frequency of the observed grating image.

φ_{O} (x, y)

represents phase which is modulated by the height of the object, and

φ_{I} (x, y)

represents the phase without the object on the reference plane [7].

Fig. 2 The optical geometry of fringe projection profilometry.

Download Full Size | PDF

By computing the 1-D Fourier transform in the x-direction of Eq. (1), the Fourier spectra of every line in the object diagram could be obtained, as is shown in Fig. 3. With a suitable filter function used, the fundamental component of the spectra could be obtained. If we use the improved Fourier transform method and acquire another image with π phase shift [19]. We could form that:

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

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

Therefore, the differences between Eq. (3) and Eq. (4) is:

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

Same processes are carried out to

I (x, y)

and the same equation can be obtained:

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

It can be seen that there is no zero component in Eq. (5) and (6) after the Fourier transform. Thus, the fundamental component can be extended toward lower frequencies near 0 and toward higher frequencies, of at least 2 times the fundamental frequency. Therefore, improved Fourier transform profilometry gives us an expanded measurement range compared to normal Fourier transform profilometry. Detailed explanations can be found in [19].

Fig. 3 Spatial frequency spectra of deformed grating pattern.

Download Full Size | PDF

Here, the Fourier Transform is applied to both Eq. (5) and (6) to obtain the fundamental component with a window function. After that, we apply the inverse Fourier transform to the fundamental component and get a complex signal [19]:

\tilde{O_{0}} (x, y) = 2 b (x, y) \exp [i 2 π f_{0} x + i φ_{O} (x, y)]

\tilde{I_{0}} (x, y) = 2 b (x, y) \exp [i 2 π f_{0} x + i φ_{I} (x, y)]

Then, the phase modulated by the height of the object can be expressed as:

Δ φ (x, y) = φ_{O} (x, y) - φ_{I} (x, y) = arc \tan \frac{Im [\tilde{O_{0}} (x, y) \tilde{I_{0}^{*}} (x, y)]}{Re [\tilde{O_{0}} (x, y) \tilde{I_{0}^{*}} (x, y)]}

According to Fig. 2, we can write that:

Δ φ (x, y) = 2 π f_{0} \overset{— —}{C D}

, so the height of the object is:

H (x, y) = \frac{l_{0} Δ φ (x, y)}{[Δ φ (x, y) - 2 π f_{0} d]}

In which

l_{0}

is the distance between the reference plane and the optical center of the CCD camera,

d

is the distance between the optical center of the CCD camera and the projector. In the real experiment,

d ≫ \overset{— —}{C D}

, so that Eq. (10) could be simplified to:

H (x, y) = \frac{l_{0} Δ φ (x, y)}{2 π f_{0} d}

Consequently, the relationship between the real height of the object and the unwrapped phase is obtained.

2.2 Peak searching algorithm for calculating the height of object

Considering the periodicity of a sinusoid, we recognize that the modulus term can be taken into account. Therefore, simply adding the sinusoids together will give a peak signal at the height of the object. We then form the following function:

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

In which,

K

is the number of the total spatial frequencies used in the measurement process,

C_{k}

is the relationship between

Δ φ (x, y)

and the real height of the object, and

h

is the variable used to find the real height of the object. According to Eq. (11),

C_{k}

here means

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

, in which

f_{k}

is the spatial frequency of the fringe pattern. Then, we have:

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

If we use only unity valued cosines and measured phase information to find the real height of the object, the actual position of the peak of the sum of cosines is system phase dependent. By using complex notation and taking the magnitude, the highly oscillating term can be removed and only the envelope remains. The envelope peak is phase independent [18]. Thus, we rewrite the function:

S (h) = \frac{1}{K} | \sum_{k = 1}^{K} \exp {[\frac{2 π f_{k} d}{l_{0}} h - Δ φ (x, y)] \cdot i} |

According to Eq. (11), we can change Eq. (14) to that:

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

In which

H

is the real height of the object. Obviously when h equals to the real height of the object, Eq. (15) exports a peak higher than any other

h

. Thus, there is a highest peak in Eq. (14) that represents the real height of the object.

As was discussed above, we have established a technique of multi-frequency fringe projection to get the height, in which the height of each point of the object is calculated independently and phase unwrapping procedure is avoided, enabling us to measure objects with large depth discontinuities.

2.3 Mathematical proof of the above methods

Considering the expression inside Eq. (15), we can write that:

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

The spatial frequencies of the fringe patterns used meet the following relationship:

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

Here,

m

is a constant. In this way, Eq. (16) can be regarded as sums of a geometric sequence. According to the summation formula of geometric series, Eq. (16) can be written as:

g (h) = \frac{\exp [\frac{2 π f_{1} d}{l_{0}} (h - H) \cdot i] \cdot {1 - \exp [\frac{2 π K m d}{l_{0}} (h - H) \cdot i]}}{1 - \exp [\frac{2 π m d}{l_{0}} (h - H) \cdot i]}

Then, applying Euler’s formula to Eq. (18), we can write that:

\begin{array}{l} g (h) = \exp [\frac{2 π f_{1} d}{l_{0}} (h - H) \cdot i] \times \frac{\exp [\frac{2 π d}{l_{0}} \times \frac{K m}{2} (h - H) \cdot i]}{\exp [\frac{2 π d}{l_{0}} \times \frac{m}{2} (h - H) \cdot i]} \\ \times \frac{\exp [\frac{2 π d}{l_{0}} \times (- \frac{K m}{2}) (h - H) \cdot i] - \exp [\frac{2 π d}{l_{0}} \times \frac{K m}{2} (h - H) \cdot i]}{\exp [\frac{2 π d}{l_{0}} \times (- \frac{m}{2}) (h - H) \cdot i] - \exp [\frac{2 π d}{l_{0}} \times \frac{m}{2} (h - H) \cdot i]} \end{array}

In order to simplify the expression above, we assume that:

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

Therefore, Eq. (19) can be written as:

g (h) = \frac{\exp (i K α)}{\exp (i α)} \times \frac{\sin (K α)}{\sin (α)} \times \exp (\frac{2 α f_{1}}{m} \cdot i) = \frac{\sin (K α)}{\sin (α)} \times \exp [(K - 1) α + \frac{2 α f_{1}}{m}] \cdot i

According to Eq. (15), Eq. (16) and the equations above, the function

S (h)

could be rewritten as:

S (h) = \frac{1}{K} | g (h) | = \frac{1}{K} \sqrt{g (h) \times g^{*} (h)} = \frac{1}{K} \times \frac{\sin (K α)}{\sin (α)}

As we know, functions that are similar to

\sin (K α) / \sin (α)

have the property that the main maximum of the function appears when

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

. This means that when

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

,

S (h)

will export the main maximum which is larger than any other positions. Equation (14) is simulated and shown in Fig. 4.

Fig. 4 The simulated curve of Eq. (14) by Matlab.

Download Full Size | PDF

We can give the following definition according to Fig. 4:

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

Resolution = \frac{l_{0}}{d K m} = \frac{l_{0}}{d (f_{k} - f_{1})}

Equation (23) is the length between the two nearest peaks and it is named ambiguity free height. Within this range, the height of the object can be found easily and accurately without encountering extra peaks. The finer the spatial frequency spacing, the greater the ambiguous free height range. Equation (24) is the full width half maximum (FWHM) of the peak and the FWHM is the nominal resolution of this method, which is related to the bandwidth of the selected frequencies as shown in Fig. 4. Therefore, the greater the bandwidth, the finer the resolution.

3. Experimental results

In the experiment, we used a DMD projector to project the fringe patterns generated by LabVIEW, and a JAI GigE CCD camera (CM-030GE) was used to collect the deformed fringe patterns. Here, we only needed to collect two images with a phase-shift of π for each frequency. The frequencies of the fringe patterns were equally spaced from 7/6 fringes/cm to 1.5 fringes/cm with a common difference of 1/30 fringes/cm.

The object used in the experiment is shown in Fig. 5. The height of the object was measured by Vernier calipers. The higher part of the object is 14.72mm high, while the lower part is 5.14mm high. The image without background intensity and higher spectrum is shown in Fig. 6(a). Figure 6(b) is the wrapped phase map of the object calculated by Fourier Transform.

Fig. 5 The object used in the experiment.

Download Full Size | PDF

Fig. 6 (a) Image without background intensity and higher spectrum, (b) Wrapped phase map calculated by Fourier Transform.

Download Full Size | PDF

After carrying out the proposed method mentioned in the last section, we show the curve of Eq. (14) for one point on the object in Fig. 7. Since we used 11 frequencies, there are 10 small peak value between the main maximums in Fig. 7.

Fig. 7 The curve of Eq. (14) for point (120,120) in the experiment.

Download Full Size | PDF

The reconstructed 3D shape of the object is shown in Fig. 8. The measured height of the higher part is 14.58mm high while the real height is 14.72mm and the lower part is 4.99mm high while the real height is 5.14mm. Therefore, for the higher part, the measurement error is 0.14mm and 0.15mm for the lower part. The accuracy of measurement is related to several factors such as the gamma effect of the projector and the CCD camera, the selected frequencies, and the resolution of the CCD camera. The accuracy could be improved with more frequencies or by improving the gray level resolution and the spatial resolution of the CCD camera.

Fig. 8 The reconstructed 3D shape of the object.

Download Full Size | PDF

The reconstruction time of the 3D shape is short, because only two images with π phase-shift are needed for each frequency. It takes about 60s to complete the measurement process.

4. Conclusion

A multi-frequency fringe projection profilometry with Fourier transform is proposed. A peak searching algorithm and its mathematical proof are presented. The measurement results of an object with large discontinuities confirm the validity of the proposed method. Moreover, the measurement range is expanded because of using Fourier Transform profilometry with π phase-shift. The proposed method requires no phase unwrapping process and each point of the object is calculated independently. Therefore, our technique is capable of measuring objects with large depth discontinuities, making it promising in industrial applications.

Acknowledgments

This work is supported by the National Natural Science Fund Committee and the Civil Aviation Administration of China Jointly Funded Project (U1233116) and Research Fund for the Doctoral Program of Higher Education of China (20120201110032)

References and links

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

2. E. Li, X. Peng, J. Xi, J. Chicharo, J. Yao, and D. Zhang, “Multi-frequency and multiple phase-shift sinusoidal fringe projection for 3D profilometry,” Opt. Express 13(5), 1561–1569 (2005). [CrossRef] [PubMed]

3. S. Vrinivasan, H. C. Liu, and M. Halioua, “Automated phase-measuring profilometry of 3-D object shapes,” Appl. Opt. 23(18), 3105–3108 (1983). [PubMed]

4. X.-Y. Su, W.-S. Zhou, G. von Bally, and D. Vukicevic, “Automated phase-measuring profilometry using defocused projection of a Ronchi grating,” Opt. Commun. 94(6), 561–573 (1992). [CrossRef]

5. H. Zhang, M. J. Lalor, and D. R. Burton, “Spatiotemporal phase unwrapping for the measurement of discontinuous objects in dynamic fringe-projection phase-shifting profilometry,” Appl. Opt. 38(16), 3534–3541 (1999). [CrossRef] [PubMed]

6. M. Takeda and K. Mutoh, “Fourier transform profilometry for the automatic measurement of 3-D object shapes,” Appl. Opt. 22(24), 3977–3982 (1983). [CrossRef] [PubMed]

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

8. 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]

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

10. L. D. Xiong and S. H. 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]

11. J. M. Huntley and H. O. Saldner, “Temporal phase-unwrapping algorithm for automated interferogram analysis,” Appl. Opt. 32(17), 3047–3052 (1993). [CrossRef] [PubMed]

12. H. O. Saldner and J. M. Huntley, “Temporal phase unwrapping: application to surface profiling of discontinuous objects,” Appl. Opt. 36(13), 2770–2775 (1997). [CrossRef] [PubMed]

13. J. M. Huntley and H. O. Saldner, “Error-reduction methods for shape measurement by temporal phase unwrapping,” J. Opt. Soc. Am. A 14(12), 3188–3196 (1997). [CrossRef]

14. J. L. Li, H. J. Su, and X. Y. Su, “Two-frequency grating used in phase-measuring profilometry,” Appl. Opt. 36(1), 277–280 (1997). [CrossRef] [PubMed]

15. Y. Wang and S. Zhang, “Superfast multifrequency phase-shifting technique with optimal pulse width modulation,” Opt. Express 19(6), 5149–5155 (2011). [CrossRef] [PubMed]

16. E. H. Kim, J. Hahn, H. Kim, and B. Lee, “Profilometry without phase unwrapping using multi-frequency and four-step phase-shift sinusoidal fringe projection,” Opt. Express 17(10), 7818–7830 (2009). [CrossRef] [PubMed]

17. C. Reich, R. Ritter, and J. Thesing, “White light heterodyne principle for 3D-measurement,” Proc. SPIE 3100, 236–244 (1997). [CrossRef]

18. 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]

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

Recovery of absolute height from wrapped phase maps for fringe projection profilometry

Abstract

1. Introduction

2. Principle and method of multi-frequency fringe projection profilometry

2.1 Obtain the wrapped phase map by Fourier transform

2.2 Peak searching algorithm for calculating the height of object

2.3 Mathematical proof of the above methods

3. Experimental results

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (8)

Equations (24)

Optics Express