Optics and Lasers in Engineering 35 (2001) 263–284
Fourier transform profilometry: a review Xianyu Su*, Wenjing Chen Opto-Electronics Department, Sichuan University, Chengdu 610064, China Received 1 September 2000; received in revised form 20 February 2001; accepted 28 February 2001
Abstract Fourier transform profilometry is one of the popular non-contact 3-D measurement methods, where a Ronchi grating or sinusoidal grating is projected onto a diffuse threedimensional surface, and the resulting deformed grating image is detected by a CCD camera and processed by a computer. This method requires only one frame (or two frames) of the deformed fringe pattern in some algorithms to retrieve the surface of measured object, so it has obvious advantage for real time data acquisition and 3-D measurement of dynamic process. In this paper, we review some algorithms in FTP, discuss some important problems, including frequency spectra overlapping, phase unwrapping, sampling, and 3-D measurement of dynamic process. With the development of computer hardware and software and availability of high-resolution image grabber, FTP method will be a promising one for acquiring 3-D data of object, and more and more researchers pay attention to it. # 2001 Elsevier Science Ltd. All rights reserved.
1. Introduction Optical 3-D non-contract profilometry has been widely used for 3-D sensing, mechanical engineering, machine vision, intelligent robots control, industry monitoring, biomedicine, dressmaking, etc. Several 3-D object profilometry methods that use structured light pattern, including Moir!e technique (MT) [1,2], phasemeasuring profilometry (PMP) [3–11], Fourier transformation profilometry (FTP) [12,13], modulation measurement profilometry (MMP) [14,15], spatial phase detection (SPD) [16–18], laser triangulation measurement [19–22], color-coded
*Corresponding author. E-mail address:
[email protected] (X. Su). 0143-8166/01/$ - see front matter # 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 1 4 3 - 8 1 6 6 ( 0 1 ) 0 0 0 2 3 - 9
264
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
fringe projection [23], gray-coded binary fringe sequences [24] have been exhaustively studied. Among them, Fourier transform profilometry, by Takeda et al. is a popular one, because of following merits, only one (or two) fringe(s) needed, full field analysis, and high precision, etc. In FTP, a Ronchi grating or a sinusoidal grating is projected onto the object surface, the depth information of the object is encoded into a deformed fringe pattern recorded by an image acquisition sensor. The surface shape can be decoded by calculating Fourier transformation, filtering in spatial frequency domain, and calculating inverse Fourier transformation. Compared with the Moir!e topography (MT), FTP can accomplish fully automatic distinction between a depression and an evaluation of the object shape. It requires no fringe order assignments or fringe center determination, and it requires no interpolation between fringes because it gives height distribution at every pixel over the entire fields. Compared with the phase-measuring profilometry (PMP) and modulation measurement profilometry (MMP), FTP requires only one or two images of the deformed fringe pattern, which makes real-time data processing and dynamic data processing possible. Whereas, PMP and MMP require many images with fixed phase variation between neighboring two images to retrieve the height distribution, and they take much time to finish phase-shifting procedure using mechanical device. So it is impossible to use them to measure dynamic object. After Takeda et al. the FTP method has been extensively studied [25–41]. A grating p phase shifting technique [28,29] can extend the measurable slope of height variation to nearly three times that of the unimproved FTP. Two-dimensional Fourier transform and 2-D hanning filtering are applied to provide a better separation of the height information from noise when speckle-like structures and discontinuities exist in the fringe pattern [31]. FTP based on TDI camera can be used to measure 3608 shape [33]. The frequency-multiplex technology [35,36] permits the 3-D shape measurement of objects that have discontinues height step and/or spatially isolated surfaces. The phase error caused by sampling in FTP is discussed in detail [40]. This paper reviews some of the developments in Fourier transform profilometry over the past years, including frequency spectra overlapping, phase unwrapping, sampling, and 3-D measurement of dynamic process, complex object phase unwrapping, 3-D phase unwrapping, etc. With the availability of high resolving CCD and high frame rate grabber, FTP has become an effective method for 3-D shape measurement.
2. Principal FTP was introduced by Takeda et al. The general geometry is shown in Fig. 1, in which the optical axes P1 P2 of a projector lens crosses that of camera lens I1 I2 at point O on a reference plane R, which is a fictitious plane normal to I1 I2 and serves as a reference, from which object height hðx; yÞ is measured. D expresses a tested point, A, C express points on R, d is the distance between P2 and I2, and L0 is the distance between I2 and O. Ronchi grating G has its lines normal to the plane of the figure, and its image is projected onto the object surface.
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
265
Fig. 1. Optical geometry.
When the measured object is put on R, the deformed fringe pattern observed through a CCD can be expressed by 1 X gðx; yÞ ¼ rðx; yÞ An expðið2pnf0 x þ njðx; yÞÞÞ: ð1Þ n¼1
When hðx; yÞ ¼ 0, the deformed grating image is written as 1 X An expðið2pnf0 x þ nj0 ðx; yÞÞÞ; g0 ðx; yÞ ¼ r0 ðx; yÞ
ð2Þ
n¼1
where rðx; yÞ and r0 ðx; yÞ are non-uniform distributions of reflectivity on the diffuse object and on the reference plane, respectively. An is the weighting factors of Fourier series, f0 is the fundamental frequency of the observed grating image, jðx; yÞ and j0 ðx; yÞ are the phase modulations resulting from the object height distribution and original phase modulation (hðx; yÞ ¼ 0), respectively. Compute the 1-D Fourier transform of Eq. (1) and obtain the Fourier spectra, as shown in Fig. 2. With a suitable filter function, the spectra are filtered to let only the fundamental component. The inverse Fourier transform is applied to the fundamental component, we obtain a complex signal g# ðx; yÞ ¼ A1 rðx; yÞ expði2pf0 x þ fðx; yÞÞ:
ð3Þ
The same operation is applied to Eq. (2), we obtain g# 0 ðx; yÞ ¼ A1 r0 ðx; yÞ expði2pf0 x þ f0 ðx; yÞÞ:
ð4Þ
The core variable varied directly with the height distribution is the phase change Djðx; yÞ: % Djðx; yÞ ¼ jðx; yÞ j0 ðx; yÞ ¼ 2pf0 AC: ð5Þ Djðx; yÞ can be obtained by Eqs. (3) and (4) as Dfðx; yÞ ¼ Imflog½#gðx; yÞ g# 0* ðx; y Þg:
ð6Þ
266
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
Fig. 2. 1-D spatial frequency spectra of deformed grating pattern.
The phase calculation by computer gives principal values ranging from p to p, and has discontinuities with 2p phase jumps. By phase unwrapping algorithm, we can obtain continuous phase distribution. Assuming the measuring system is perfect and the system parameters (L0 and d) have no error, then height distribution can be computed by the following formula: L0 Dfðx; yÞ : ð7Þ hðx; yÞ ¼ Dfðx; yÞ 2pf0 d Substituting Eq. (6) into Eq. (7), hðx; yÞ is obtained. When the measuring system is not perfect, how to solve the problems have been discussed [21,25]. Since Fourier spectra of the deformed fringe pattern usually have multiple components, which can be expressed as n @jðx; yÞ fn ¼ nf0 þ : ð8Þ 2p @x In order to restore the object surface correctly, the fundamental component must be separated from all other spectra, as shown in Fig. 2, it is necessary that ðf1 Þmax 5ðfn Þmin ;
n > 1;
ðf1 Þmin > fb :
ð9Þ ð10Þ
Consequently, the phase variation caused by height modulation must be limited in @fðx; yÞ ð11Þ @x 52pf0 =3: max Assuming that L0 4hðx; yÞ, we have approximately 2pf0 d hðx; yÞ: fðx; yÞ Dfðx; yÞ ¼ L0 Substituting Eq. (12) into Eq. (11), we obtain @hðx; yÞ @x 5L0 =3d: max
ð12Þ
ð13Þ
That is: FTP can be used only for surfaces on which the slopes do not exceed this limitation. When the measurable slope of height variation extend the limitation, the fundamental components overlap the zero component and other high components.
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
267
Aliasing will influence the precision of FTP, even correct reconstruction cannot be obtained. In following section, we review the techniques used to reduce frequency overlapping and improve the measurement precision.
3. Several improvements in FTP 3.1. Quasi-sine projection and p phase shifting technique [28,29] In FTP, we apply quasi-sine projection technique and p phase shifting technique to merely make the fundamental component exist in spatial frequency domain. In this case, the lower frequency part of the fundamental component can extend to zero, and the higher part can extend to 2f0 without overlaps. FTP can measure the object with higher variation. The intensity distribution on the object surface can be expressed as gðx; yÞ ¼ aðx; yÞ þ bðx; yÞ cosð2pf0 x þ fðx; yÞÞ:
ð14Þ
To eliminate the influence of background intensity aðx; yÞ, the optical field is sampled twice. At the second sampling, the grating is moved half of the grating period, and the other conditions remain unchanged. The difference between the two fringes is gðx; yÞ ¼ 2bðx; yÞ cosð2pf0 x þ fðx; yÞÞ
ð15Þ
where no zero component and higher components exist in frequency domain. So the fundamental component can be extended toward lower frequency, nearly to 0, and toward higher frequency, at least to nearly 2f0 in the spectral domain. The results of quasi-sine projection and p phase shifting technique result in the larger range of the measurement, that is @fðx; yÞ ð16Þ @x 52pf0 ; max @hðx; yÞ ð17Þ @x 5L0 =d: max It is obvious that the improved method can raise the measurable slope of the height variation nearly three times while other system parameters remain unchanged. A primary experiment is used to verify our demonstration. The measured object is a triangle, which base is 40 mm, and the height is 60 mm, the slopes of the triangle are 3 and –3. The parameter of the system is L0 =d ¼ 5:2. Since L0 =3d5@hðx; yÞ=@xmax 5L0 =d, Fig. 3(a) shows the frequency spectra of the deformed fringe, the overlap is obvious. So the two sides of the triangle cannot be restored completely by unimproved FTP, as shown in Fig. 3(b). Using quasi-sine projection and p phase shifting technique, the background intensity aðx; yÞ and higher spectra are eliminated, as shown in Fig. 3(c), the surface of the triangle can be completely restored, as shown in Fig. 3(d).
268
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
Fig. 3. (a) Spatial frequency spectra of deformed grating pattern and restored shape by unimproved FTP. (b) Spatial frequency spectra of deformed grating pattern and restored shape by improved FTP.
3.2. Two-dimensional Fourier transform profilometry [31] To accomplish automatic measurement of a coarse object, on which there are discontinuities and speckle-like structures in the fringe pattern, for example, sandbody, two-dimensional Fourier transform profilometry is applied. In this case, the deformed fringe pattern with noise can be expressed as gðx; yÞ ¼ rðx; yÞ
1 X
bn expðið2pnf0 x þ njðx; yÞÞÞ þ Nðx; yÞ;
ð18Þ
n¼1
where Nðx; yÞ represents 2-D noise distribution. The 2-D Fourier transform is applied to compute the spectra of Eq. (18), then 2-D hanning filter function, which is expressed by Eq. (19), is used to reinforce the frequencies around the carrier frequency f0 and attenuate the rest more as the distance from f0 is increased. The frequencies caused by speckle-like structure and the discontinuity are eliminated,
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
269
Fig. 3. (Continued.)
finally 2-D inverse Fourier transform is applied to the fundamental component. fy fx f0 Hðfx ; fy Þ ¼ 0:25 1 þ cos bp 1 þ cos bp : ð19Þ fcx fyx The same operation is applied to the reference fringe. Applying the phase algorithm and phase unwrapping algorithm, we obtain the 2-D continuous phase distribution. Fig. 4(a) shows the surface of a Chinese word ‘‘jiang’’ written on a sand-body, on which a Ronchi grating pattern is projected. By 2-D FTP method, the shape of the object can be restored correctly, as shown in Fig. 4(b). 3.3. Fourier transform profilometry based on TDI camera [33] Time delay and integration (TDI) camera can be easily adapted for recording image in a high speed. Combined with TDI camera, Fourier transform profilometry
270
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
Fig. 4. (a) A Chinese word ‘‘jiang’’ shape on the sand with a grating pattern projected onto it. (b) The restored surface of ‘‘jiang’’ using 2-D FTP.
can be used for 3608 shape measurement. When the object rotates for 3608, TDI camera records single stripe structured light recorded sequentially on one image. The final image contains the surface information of the 3608 object. The fringe pattern is processed by 2-D Fourier transform profilometry to obtain the phase change caused
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
271
Fig. 5. The schematic for the 3608 profilometry.
by height distribution of the 3608 object. Then it is easy to calculate the actual height distribution. Since the sampling and processing are separated, and only one image is needed, the speed of the method can be improved rapidly. The schematic of the 3608 shape testing is given in Fig. 5. Illumination is by a thin line of light from a laser diode. The TDI camera is mounted so that the scanning direction is perpendicular to the rotation axis of the test object. When the object has been rotated for 3608, a distorted fringe pattern is formed on a single image. As a demonstration, a mannequin head is used. The fringe pattern image of a 3608 light stripe obtained by TDI camera is shown in Fig. 6(a). Fig. 6(b) is the contour map with gray scale presentation. Fig. 6(c) is the reconstructed profile of the mannequin head. 3.4. Frequency-multiplex Fourier transform profilometry [35,36] In 1997, Takeda et al. proposed frequency-multiplex Fourier transform profilometry based on spatial frequency multiplexing combined with the Gushov– Solodkin phase unwrapping algorithm. The technique permits the three-dimensional shape measurement of objects that have discontinuous height steps and/or spatially isolated surface, which has not possible by conventional FTP. An important feature of the technique is that it requires only a single fringe pattern. So it is suitable for instantaneous 3-D shape measurement of discontinuous objects in fast motion. 4. The influence of sampling in FTP [40] Employing continuous Fourier transform analytical method, Takeda et al. analyzed the problems that existed in FTP, such as frequency aliasing, measurable range, etc. [13]. But in practice, the deformed fringe pattern captured by CCD
272
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
Fig. 6. (a) Projected grating lines on unwrapped mannequin head. (b) The contour map with gray-scale presentation. (c) 3-D display of the mannequin.
camera is discrete, and the discrete Fourier transform (DFT) is applied. The obvious shortcoming in theory analysis results in the inconsistency between theory and experience, especially if high-order spectra existed. To obtain the correct reconstruction of the measured object, the influence of sampling in FTP must be discussed. The deformed fringe pattern caused by CCD camera can be expressed as x y gðx; yÞ ¼ qn ðx; yÞ expði2pnf0 xÞcomb comb ; Dx Dy n¼1 1 X
ð20Þ
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
273
where qn ðx; yÞ ¼ rðx; yÞbn expðinjðx; yÞÞ:
ð21Þ
Dx and Dy denote sampling intervals in x, y directions. For simplicity, we compute the 1-D Fourier transform of Eq. (20) for the variable x only, with y fixed, and obtain 1 1 X X y Gs ðf ; yÞ ¼ Qn ðf nf0 Nmf0 ; yÞ comb ; ð22Þ Dy N¼1 n¼1 P where Nmf0 ¼ NDf ¼ N=Dx, Eq. (22) states that the function 1 n¼1 Qn ðf nf0 ; yÞ appears in replica at intervals of mf0 in both directions along x-axis, which is called a general spectrum region. To eliminate the overlapping between spectral regions, m must beP big enough. As shown in Fig. 7, the dark lines correspond to the right-hand side of 1 nf0 ; yÞ, the light lines correspond to the left-hand side of its n¼1 Qn ðf P 1 adjacent frequency n¼1 Qn ðf nf0 mf0 ; yÞ. According to the definition of a local spatial frequency fn for the nth spectrum component, shown as Eq. (8), and assuming the fringe pattern to be a band-limited signal, then ðfn Þmax corresponds to the cutoff frequency of the Fourier spectra. Only whenPðfn Þmax does not exceed 1 mf0 ðfn Þmax , does the frequency P1 spectra determined by n¼1 Qn ðf nf0 ; yÞ not overlap those determined by n¼1 Qn ðf nf0 mf0 ; yÞ. In other words, ðfn Þmax 5 mf0 ðfn Þmax must be satisfied. In a word, in FTP, the following three inequations must be satisfied: ðf1 Þmax 5ðfn Þmin ;
n > 1;
ð23Þ
ðf1 Þmin > fb ;
ð24Þ
ðf1 Þmax 5mf0 ðfn Þmax :
ð25Þ
Substituting Eq. (8) into the above inequations, we obtain @hðx; yÞ L0 @x 53d ; max
Fig. 7. Replicated spectrum island.
ð26Þ
274
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
@hðx; yÞ m n 1 L0 @x 5 n þ 1 d ; max
ð27Þ
where m ¼ Df =f0 , and the bigger m results in the larger range of the measurement. Since the weighing factor of the second-order frequency spectrum is much bigger than that of other spectra, the phase error caused by overlapping between the fundamental frequency and the second-order spectrum in the same island is larger than that from the adjacent period. That is, only when Eq. (26) is satisfied, does increasing sampling frequency to decrease overlapping between periods have practical meaning for decreasing the phase errors. It is necessary that 1 L0 m n 1 L0 5 ; 3 d nþ1 d
ð28Þ
that is m543ðn þ 1Þ:
ð29Þ
So owing to the discrete nature of the DFT, in practical measurement, we must ensure that the fundamental component does not overlap the zero component and other higher-order spectra in the same Frequency Island. Furthermore, we must also ensure that the fundamental component does not overlap the higher spectra from the adjacent frequency islands. is used in our its base line is 54 mm and height is 27 mm, A triangle experiment, @h=@x ¼ 1, L0 =d ¼ 4:5, @h=@x 5L0 =3d. According to the theories proposed by max previous authors [12,13], the triangle surface can be correctly restored by FTP. But taken into account the overlapping between adjacent frequency period caused by the third spectrum of Ronchi grating, m > 16 3 5:3, according to Eq. (29). If the sampling frequency is 4f0 , correct reconstruction cannot be obtained. The frequency spectra and the restored surface are shown in Fig. 8(a). If the sampling frequency is 7f0 , we can obtain the correct reconstruction. The frequency spectra and the restored surface are shown in Fig. 8(b).
5. Dynamic 3-D shape measurement method based on FTP Now we review the method based on Fourier transform profilometry for dynamic 3-D shape measurement. First, we record a fringe pattern for hðx; yÞ ¼ 0: g0 ðx; yÞ ¼
1 X
An expðið2pnf0 x þ nj0 ðx; yÞÞÞ;
ð30Þ
n¼1
then put a dynamic 3-D object in the optical field, and record and store a sequence of deformed fringe patterns rapidly. The intensity distributions of these fringe patterns
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
275
can be expressed as 1 X
gðx; y; tÞ ¼
An rðx; y; tÞ expðið2pnf0 x þ njðx; y; tÞÞÞ
ðt ¼ 1; 2; . . . ; sÞ
n¼1
ð31Þ in Eqs. (30) and (31), rðx; y; tÞ and jðx; y; tÞ represent a non-uniform distribution of reflectivity and phase modulation resulting from the object height distribution at the time of t, respectively, which are a slowly varying function relating to position and time. s denotes the total frame number of the patterns recorded by CCD. To obtain jðx; y; tÞ separately from the unwanted amplitude variation rðx; y; tÞ caused by non-uniform reflectively on object surface. We compute the Fourier transform of deformed fringes only for the variable x, with y being fixed Gðf ; y; tÞ ¼
1 X
Qn ðf nf0 ; y; tÞ;
ð32Þ
n¼1
where Gðf ; y; tÞ and Qn ðf ; y; tÞ are the 1-D Fourier spectra of gðx; y; tÞ and qn ðx; y; tÞ, qn ðx; y; tÞ ¼ bn rðx; y; tÞ expðinfðx; y; tÞÞ, respectively. If the fundamental component separates from the other spectra, using a suitable filter function the spectra can be filtered to let only fundamental component through. The inverse Fourier transform is applied to the fundamental component, we get g# ðx; y; tÞ ¼ A1 rðx; y; tÞ expði2pf0 x þ fðx; y; tÞÞ:
ð33Þ
The same operation is done on the reference fringe to obtain g# 0 ðx; yÞ ¼ A1 r0 ðx; yÞ expði2pf0 x þ f0 ðx; yÞÞ:
ð34Þ
There are two methods to obtain the dynamic phase change resulting from height distribution. 5.1. Direct phase algorithm Calculating the multiplication of the conjugation of the g# 0 ðx; yÞ with g# ðx; y; tÞ, we obtain g# ðx; y; tÞ#g0* ðx; yÞ ¼ jA1 j2 rðx; y; tÞ expðiDfðx; y; tÞÞ:
ð35Þ
The phase change Djðx; y; tÞ can be obtained by Djðx; y; tÞ ¼ jðx; y; tÞ j0 ðx; yÞ ¼ arctg
Im½#gðx; y; tÞ#g0* ðx; yÞ ; Re½#gðx; y; tÞ#g0* ðx; yÞ
ð36Þ
where Im[ ] and Re[ ] represent the imaginary part and real part, respectively. By this method, we can obtain a sequence of phase maps, which include the fluctuation information of dynamic object.
276
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
5.2. Phase difference algorithm By calculating the product of the g# 0* ðx; yÞ with g# ðx; y; 1Þ and that of g# * ðx; y; t 1Þ with g# ðx; y; tÞ in neighboring time, we obtain g# ðx; y; 1Þ#g0* ðx; yÞ ¼ jA1 j2 rðx; y; 1Þ expðiDf1 ðx; yÞÞ;
ð37Þ
g# ðx; y; tÞ#g * ðx; y; t 1Þ ¼ jA1 j2 rðx; y; tÞrðx; y; t 1Þ expðiDft ðx; yÞÞ:
ð38Þ
Dft ðx; yÞ can be expressed as 8
Im g# ðx; > > > fðx; y; tÞ f0 ðx; yÞ ¼ arctg
< Re g# ðx; Dft ðx; yÞ ¼ > Im½g# ðx; > > : fðx; y; tÞ f0 ðx; y; t 1Þ ¼ arctg Re½g# ðx;
y; y; y; y;
1Þ#g0* ðx; yÞ 1Þ#g0* ðx; yÞ tÞ#g * ðx; y; t 1Þ : tÞ#g * ðx; y; t 1Þ
ðt ¼ 1Þ ðt > 1Þ
ð39Þ By computing the finite sum of Eq. (39), we get t X Dfn ðx; yÞ ¼ jðx; y; tÞ j0 ðx; yÞ:
ð40Þ
n¼1
Comparing Eq. (36) to Eq. (40), we get t X Dfn ðx; yÞ: Dfðx; y; tÞ ¼
ð41Þ
n¼1
The phase distribution in any time computed by the two methods discussed above lies in (p p). To obtain continuous phase distribution, we must unwrap the phase along x, y coordinates and time dimension t. In evidence, the phase difference of fðx; y; tÞ fðx; y; t 1Þ is always smaller than that of fðx; y; tÞ f0 ðx; yÞ, this provides a new approach for phase unwrapping. 6. Phase unwrapping in FTP method Phase unwrapping is a critical but challenging step in any grating projection profilometry, including Fourier transform profilometry. In recent years, phase unwrapping problem has been widely studied [42–56]. In this paper, we will concentrate on the phase unwrapping based on the modulation ordering and extend 2-D phase unwrapping algorithm into 3-D phase space for dynamic 3-D measurement. 6.1. Phase unwrapping based on modulation for complex object [5,45–51] Since the phase calculation by any inverse trigonometric function provides the principal values ranging from p to p, for phase data without noise, these discontinuities can be corrected easily by adding or subtracting 2p according to the phase jump ranging from p to p or vice versa. This is called the phase unwrapping
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
277
Fig. 8. Measured object’s spectra and its retrieved shape. (a) Measured object’s spectra and it’s retrieved shape when m=4. (b) Measured object’s spectra and its retrieved shape when m=7.
procedure. If the wrapped phase is reliable everywhere and the maximum phase variation between neighboring pixels is less than p everywhere, the phase unwrapping procedure consists of adding multiples of 2p to the phase when discontinuities appear. The unwrapped phase is correctly written as fu ðx; yÞ ¼ fðx; yÞ þ 2kp:
ð42Þ
But in practice, phase unwrapping is a difficult problem for complex shape measurement. Many factors, such as noise, local shadows, under-sampling, fringe discontinuities and irregular surface brightness make the actual unwrapping procedure complicated and path-dependent. Combined with modulation analysis technique, we propose a phase unwrapping algorithm suitable for complex 3-D shape measurement. First, we define modulation function as mðx; yÞ ¼ j#gðx; yÞj ¼ A1 rðx; yÞ;
ð43Þ
which means the modulation is direct ratio to non-uniform distributions of reflectivity rðx; yÞ, obviously, in the areas of local shadow and in the regions of
278
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
abrupt surface discontinuities, modulation is low, which means that the phase is uncertain. Therefore, the modulation function can be selected as one of the criteria of reliability. Phase unwrapping based on modulation has been successfully used in FTP for complex object [51]. In FTP, filtering by the digital weighting filter in frequency domain is carried out to make modulation function contain information about fringe density. For example, by hanning filter function, we can reinforce the frequency components around carrier frequency, and attenuate the parts more as the distance from f0 is increased. This associates the modulation function with fringe density. The phase unwrapping based on modulation start from the pixel with higher intensity modulation. First, we produce a pixel queue by comparing the intensity modulation of each pixel in 3 * 3 adjacent region of start point, so that the pixel with the highest modulation is put on the top of the queue. Then, the top pixel is pulled out for unwrapping and meanwhile new pixels in 3 * 3 adjacent region of the pixel are put on the queue according to the modulation ordering. When the phase unwrapped areas are more than one pixel, the ordering operation is done for all the pixels on outer layer of boundary of the phase unwrapped areas, so that the pixel with the highest modulation is put on the top of the queue. The advantage of this approach is that the path of phase unwrapping is always along the direction from the pixel with higher modulation to the pixel with lower modulation until phase is unwrapped. We have verified our algorithm by experiment. The measured object is face of a model cat with an abrupt region around the mouth, so it is difficult to obtain the correct reconstruction. Employing the phase unwrapping based on modulation yields a satisfactory result. Fig. 9(a) shows an intermediate step of phase unwrapping, the bright region represents the unwrapped part with higher modulation, the dark region denotes that the part with lower modulation will be unwrapped later. The unwrapping path successfully goes round the shadows, fringe discontinuity and possible defects. Fig. 9(b) is the result of the reconstruction of the complex 3-D object. 6.2. 3-D phase unwrapping for dynamic object 6.2.1. Direct 3-D phase unwrapping In dynamic 3-D surface measurement, the wrapped phase maps computed is a function about position (x, y) and time (t). If it is reliable everywhere and the phase difference is less than p between neighboring pixels in 3-D phase space, the unwrapping problem is trivial. The phase unwrapping can be carried out along any direction. Fig. 10 shows one of the possible unwrapping path. First, we select any pixel (i.e. point O) at initial time (t ¼ 1) as the start point to be unwrapped. 1-D phase unwrapping procedure is carried out along the t direction. Afterwards, columns from the unwrapped pixels in every phase map are selected to be unwrapped. At last, all rows from the unwrapped column in every phase map are unwrapped, and natural phases in the whole 3-D space are obtained. Of course, column and row operations
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
279
Fig. 9. (a) Intermediate stage of phase unwrapping based on modulation ordering. (b) 3-D reconstruction of complex object with gray-scale presentation.
can be interchanged. An obvious disadvantage of direct phase unwrapping is that when the wrapped phase is not reliable in a fixed time we cannot obtain correct reconstruction.
280
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
Fig. 10. Scheme of the 3-D direct phase unwrapping algorithm.
6.2.2. 3-D phase unwrapping based on modulation ordering We define modulation function as mðx; y; tÞ ¼ j#gðx; y; tÞj ¼ A1 rðx; y; tÞ;
ð44Þ
where mðx; y; tÞ represents modulation in every pixel and in difference time. Now we extend the 2-D phase unwrapping procedure based on modulation to 3-D phase space. The unwrapping procedure is same as that of the 2-D phase unwrapping except comparing the intensity modulation of each pixel in 3 * 3 * 3 adjacent region. The unwrapping scheme is as shown in Fig. 11, the gray part represents start pixel for unwrapping, lines with arrow display one of the phase unwrapping paths. 6.2.3. 3-D phase unwrapping based on phase difference algorithm Eq. (41) indicates a new 3-D phase unwrapping algorithm. The scheme is shown in Fig. 12, in which Dft1 and Dft; t1 represent the unwrapped phase by calculating phase sum from t=1 to t1 and wrapped phase difference between t and t-1, respectively. The feature of this method is that only one frame wrapped phase with reliable phase everywhere (i.e. the initial phase map corresponds to t ¼ 1) is needed to unwrap by 2-D phase unwrapping, and the unwrapping procedure is very simple in this condition. The phase values in the whole 3-D phase field can be obtained by calculating the sum along the t direction. When frame rate of frame grabber is high enough, the time interval between two frame fringe patterns grabbed is small, resulting in phase difference between neighboring pixels being less than p in the t direction. So 3-D phase unwrapping based on phase difference algorithm will play an important role in dynamic 3-D measurement. . A primary experience is applied in our demonstration. We measure the variation of human chest while breathing, as shown in Fig. 13. The process of data collection costs about 4 second. The time interval is 0.2 s, and the interval of contour plot is 1 mm. Whereas, PMP method needs to cost more than 20 s to collect the data at one time while not breathing.
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
281
Fig. 11. Scheme of the unwrapping path based on modulation ordering algorithm in 3-D space.
Fig. 12. Scheme of the unwrapping path based on phase difference algorithm in 3-D space.
About phase unwrapping, there are other algorithms, for example, additional gray-coded fringe sequence can be used to identify the fringe order to obtain the natural phase [23].
7. Conclusions As a non-contract measurement method, FTP method will play an important role in 3-D sensing field in the future because of only one or two images of the deformed fringe needed. Especially, with the availability of high-resolution CCD, for example, up to 5k * 5k pixels CCD sensor are commercially available. Combined with macroscanning technique, the resolution of CCD can be as high as 20k * 20k. So we can project high-density-grating image onto the measured object to separate the fundamental component from other unwanted spectra and capture the deformed fringe pattern by high-resolution CCD to avoid frequency aliasing between adjacent period. In future, the measurement precision of FTP can be increased largely.
282
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
Fig. 13. The contour map of human chest breathing, the time interval is 0.2 s, the interval of contour plot is 1 mm.
Since FTP method needs only one or two deformed fringe patterns to acquire the measured object surface, it has wide application in areas such as high-speed quality control and robot vision, etc. This paper reviews some of the developments in Fourier transform profilometry over the past years, including frequency spectra overlapping, phase unwrapping, sampling, and 3-D measurement of dynamic process, complex object phase unwrapping, 3-D phase unwrapping, etc. With the development of computer hardware and software and availability of high-resolution image grabber, FTP method will be a promising one for acquiring 3-D data of object, and more and more researchers pay attention to it. Acknowledgements This project was supported by the National Natural Science Foundation of China.
References [1] Takasaki H. Generation of surface contours by moir!e pattern. Appl Opt 1970;9(4):942–7. [2] Yoshizawa T. The recent trend of moir!e metrology. J Robustic Mech 1991;3(3):80–5. [3] Srinivasan V, Liu HC, Halioua M. Automated phase measuring profilometry of 3-D diffuse object. Appl Opt 1984;23(18):3l05–8. [4] Su X-Y, Zhou W-S, von Bally V, Vukicevic D. Automated phase-measuring profilometry using defocused projection of a Ronchi grating. Opt Commun 1992;94(6):561–73. [5] Su X-Y, von Bally G, Vukicevic D. Phase-stepping grating profilometry: utilization of intensity modulation analysis in complex objects evaluation. Opt Commun 1993;98(1):141–50.
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
283
[6] Zhou W-S, Su X-Y. A direct mapping algorithm for phase-measuring profilometry. J Mod Opt 1994;41(1):89–93. [7] Li JL, Su X, Su HJ. Removal of carrier frequency in phase-shifting techniques. Opt Lasers Eng 1998;30(1):107–15. [8] Li J, Su X, Su H. New PMP method with two-frequency grating projection. Proc SPIE 1996;2778:481–2. [9] Lian X, Su X. Computer simulation of a 3-D sensing system with structured illumination. Opt Laser Eng 1997;27:379–93. [10] Li J, Su H, Su X. Two-frequency grating used in phase-measuring profilometry. Appl Opt 1997;36(1):277–80. [11] Zhang H, Lalor MJ, Burton DR. Spatio temporal phase unwrapping for the measurement of discontinuous objects in dynamic fringe-projection phase-shifting profilometry. Appl Opt 1999;38(16):3534–41. [12] Takeda M, Ina H, Koboyashi S. Fourier-transform method of fringe-pattern analysis for computerbased topography and interferometry. J Opt Soc Am 1982;72(1):l56–60. [13] Takeda M, Motoh K. Fourier transform profilometry for the automatic measurement of 3-D object shapes. Appl Opt 1983;22(24):3977–82. [14] Su X, Su L. New 3D profilometry based on modulation measurement. Proc SPIE 1998;3558:1–7. [15] Su L, Su X, Li W. Application of modulation measurement profilometry to objects with surface holes. Appl Opt 1999;38(7):1153–8. [16] Toyooka S, Iwasa Y. Automatic profilometry of 3-D diffuse objects by spatial phase detection. Appl Opt 1986;25(10):30l2–8. [17] Sajan MR, Tay CJ, Shang HM, et al. TDI imaging}a tool for profilometry and automated visual inspection. Opt Lasers Eng 1998;29(6):403–11. [18] Sajan MR, Tay CJ, Shang HM, Asundi A. Improved spatial phase detection for profilometry using a TDI imager. Opt Commun 1998;150(1-6):66–70. [19] Cheng X-X, Su X-Y, Gou L-R. Automatic measurement method for 3608 profilometry of 3-D diffuse objects. Appl Opt 1991;30(10):1274–8. [20] Li J-L, Su X-Y, Zhou W-S. The 3-D sensing using laser sheet projection: influence of speckle Opt Rev 1995;2(2):144–9. [21] Asundi A, Zhou W-S. Unified calibration technique and its applications in optical triangular profilometry. Appl Opt 1999;38(16):3556–61. [22] Asundi A, Zhou W-S. Mapping algorithm for 360-deg profilometry with time delayed integration imaging. Opt Eng 1999;38(2):339–44. [23] Hausler G, Ritter D. Parallel three-dimension sensing by color-coded triangulation. Appl Opt 1993;32(35):7164–9. [24] Besi PJ, et al. Active optical image sensors. In: Sanz JLC, editor. Advances in machine vision. Berlin: Springer, 1989. [25] C’hen Y, Schenk AF. A rigorous calibration method for digital cameras. Int Arch Photogrammetry Remote Sensing 1992;27:199–205. [26] Kujawifiska M, Spik A, Wojciak J. Fringe pattern analysis using Fourier transform techniques in Interferometry ‘89: 100 Years afer Michelson. Proc SPIE 1989;1121:l30–5. [27] Kujawinska M, Wojciak J. High accuracy Fourier transform fringe pattern analysis. Opt Lasers Eng 1991;l4:325–39. [28] Su X-Y, Li J, Gou L-R. An improved Fourier transform profilometry. Proc SPIE 1988;954:32–5. [29] Li J, Su X-Y, Gou L-R. An improved Fourier transform profilometry for automatic measurement of 3-D object shapes. Opt Eng 1990;29(12):1439–44. [30] Bone DJ. Fourier fringe analysis: the two-dimensional phase unwrapping problem. Appl Opt 1991;30(25):3627–32. [31] Lin J-F, Su X-Y. Two-dimensional Fourier transform profilometry for the automatic measurement of three-dimensional object shapes. Opt Eng 1995;34(11):3297–302. [32] Su X, Su L, Li W. A new Fourier transform profilometry based on modulation measurement. ICO18, 1999.
284
X. Su, W. Chen / Optics and Lasers in Engineering 35 (2001) 263–284
[33] Su X, Sajan MR, Asundi A. Fourier transform profilometry for 360-degree shape using TDI camera. International Conference on Experimental Mechanics Advances and Applications Singapore, December 4–6,1996. [34] Yi J, Huang S. Modified Fourier transform profilometry for the measurement of 3-D steep shapes. Opt Lasers Engng 1997;27(5):493–505. [35] Takeda M, Gu Q, Kinoshita M, et al. Frequency-multiplex Fourier-transform profilometry: A singleshot three-dimensional shape measurement of objects with large height discontinuities and/or surface isolations. Appl Opt 1997;36(22):5347–54. [36] Hao YD, Zhao Y, Li DC. Multifrequency grating projection profilometry based on the nonlinear excess fraction method. Appl Opt 1999;38(19):4106–10. [37] Burton DR, Goodall AJ, Atkinson JT, et al. The use of carrier frequency-shifting for the elimination of phase discontinuities in fourier-transform profilometry. Opt Laser Eng 1995;23(4):245–57. [38] FemAndez A, Kaufmann GH, Doval AF, et al. Comparison of carrier removal methods in the analysis of TV holography fringes by the Fourier transform method. Opt Eng 1998;37(ll):2899–905. [39] Pandit SM, Chan DP. Comparison of Fourier-transform and data-dependent system profilometry by use of interferometric regeneration. Appl Opt 1999;38(19):4095–102. [40] Chen W, Yang H, Su X. Error caused by sampling in Fourier transform profilometry. Opt Eng 1999;38(6):927–31. [41] Kozloshi J, Serra G. Analysis of the complex phase error introduced by the application of Fourier transform method. Journal of Modern Optics 1999;46(6):957–71. [42] Itoh K. Analysis of phase unwrapping algorithm. Appl Opt 1982;21(14):2470. [43] Huntly JM. Noise immune phase unwrapping algorithm. Appl Opt 1989;28:3268–70. [44] Judge TR, Bryanston-Cross PJ. A review of phase unwrapping techniques in fringe analysis. Opt Lasers Eng 1994;21(4):199–239. [45] Su X-Y, Zhou W-S. Complex object profilometry and its application for dentistry. Proc SPIE 1994;2132:484–9. [46] Su X-Y, Zarubin AM, von Bally G. Modulation analysis of phase-shifted holographic interferograms. Opt Commun 1994;105(5):379–87. [47] Zhou X, Su X-Y. Effect of modulation transform function of a digital acquisition device on phasemeasuring profilometry. Appl Opt 1994;33(35):8210. [48] Li J-L, Su X-Y, Zhou W-S. The 3-D sensing using laser sheet projection: influence of speckle. Opt Rev 1995;2(2):144. [49] Su X. Phase unwrapping techniques for 3-D shape measurement (invited paper). Proc SPIE 1996;2866:460–5. [50] Li J, Su X, Li J. Phase Unwrapping algorithm-based on reliability and edge-detection. Opt Eng 1997;36(6):1685–90. [51] Su X, Xue L. Phase unwrapping algorithm based on fringe frequency analysis in Fourier transform profilometry. Opt Eng 2001;40(4). [52] Su H, Li J, Su X. Phase algorithm without the influence of carrier frequency. Opt Eng 1997;36(6):1799–805. [53] Ghiglia DC, Mastin GA, Romero LA. Cellular-automata method for phase unwrapping. J Opt Soc Am A 1987;4(l):267–80. [54] Spik A, Robinson DW. Investigation of the cellular automata method for phase unwrapping and its implementation on an array processor. Opt Lasers Eng 1991;14(l):25–37. [55] Takeda M, Abe T. Phase unwrapping by a maximum cross-amplitude spanning tree algorithm: a comparative study. Opt Eng 1996;35(8):2345–51. [56] Cuevas FJ, Servin M, Rodriguez-Vera R. Depth object recovery using radial basis functions. Opt Commun 1999;163(4–6):270–7.