Analysis and identification of phase error in phase measuring profilometry

Feng Chen; Xianyu Su; Liqun Xiang

doi:10.1364/OE.18.011300

1. Introduction

3D sensing by non-contact optical methods has been extensively studied for applications in automated manufacturing, quality control, biomedicine, machine vision, entertainment, etc. Phase measuring profilometry (PMP) [1] is an optical 3D sensing technology using fringe pattern projecting and phase-shifting. It can retrieve 3D data of object rapidly and precisely with low requirement of equipments. Research work on this subject was mainly focused on establishing effective phase-measuring algorithms and phase unwrapping algorithm [2-7].

The charge-coupled device (CCD) is used to convert the deformed fringe pattern into digital image in practical measurement. The photo-behavior of CCD greatly affects the accuracy of PMP. Eryi Hu discussed the phase-recovering algorithm for saturated fringe patterns [8,9]. In practice, the phase is quite unreliable in shadow, area with low modulation and the abrupt discontinuities because of the illumination, modulation and sampling problems. Some existing methods, such as intensity modulation, spatial frequency etc. [10], are efficient in the two former cases. They generate a reliable mask to identify the bad data which are caused by the modulation or illumination problem. But these methods do not efficiently work at the abrupt discontinuities. We analyze the effect of CCD sampling on PMP by viewing the phase modulated by the height of object is the argument of a vector and introduce temporal phase unwrapping algorithm (TPU) [11,12] into the identification of the abrupt discontinuities. The effect of CCD sampling on PMP and the identification method are presented in section 2 and 3 separately. In section 4, we present our computer simulation and experiment. The conclusion is included in the last section.

2. Effect of CCD sampling on PMP

The optical geometry of PMP system is shown in Fig. 1 . If projector and the imaging system are both telecentric optical systems, the relationship between height and phase is

Fig. 1 The optical geometry of PMP system.

Download Full Size | PDF

ϕ (x, y) = 2 π h / λ_{e}

Where $λ_{e} = P_{0} / \tan θ$ is called the equivalent wavelength. The equivalent wavelength $λ_{e}$ is determined by the spatial period $P_{0}$ of the sinusoidal grating on reference and the angle θ between principal axes of projector and that of the imaging system. The wrapped phase is obtained by using the Eq. (2) in the N-step phase calculation algorithm.

ϕ (x, y) = atan [\sum_{n = 1}^{N} I_{n} (x, y) \sin (2 n π / N) / \sum_{n = 1}^{N} I_{n} (x, y) \cos (2 n π / N)]

I_{n} (x, y)

is the intensity distribution of the deformed fringe pattern when the grating with separate phase-shifts

2 n π / N, (n = 1, \dots N)

is projected onto the object surface. In this paper,

ϕ (x, y)

is viewed as an argument of a vector

\vec{V} (x, y)

, Eq. (3).

\vec{V} (x, y) = [\sum_{n = 1}^{N} I_{n} (x, y) \cos (2 n π / N)] \vec{i} + [\sum_{n = 1}^{N} I_{n} (x, y) \sin (2 n π / N)] \vec{j}

When the deformed fringe pattern on the object is captured by CCD, the continuous 2D image is converted into a discrete digital image. The value of each pixel is in direct relation to the integral intensity of a corresponding cell in the CCD plane (Fig. 2 ). The relationship between intensity of continuous 2D image $I_{n} (x, y)$ and the value of pixel $I_{n} (i, j)$ is shown as Eq. (4) where we simply assume that $I_{n} (i, j)$ is linear to the integration of $I_{n} (x, y)$ and the factor is simplified as 1 because the quantization error is not the main issue we concern.

I_{n} (i, j) = \int \int_{x_{1 i}, y_{1 j}}^{x_{2 i}, y_{2 j}} I_{n} (x, y) d x d y

where i and j are discrete coordinates in the image.

I_{n} (i, j)

denotes the value of the pixel in ith column and jth row. The value of the pixel

(i, j)

is in direct relation to the integral intensity in rectangle

(x_{1 i} : x_{2 i}, y_{1 j}, y {}_{2 j})

. Therefore the wrapped phase in discrete form is

\begin{matrix} ϕ (i, j) = atan \frac{\sum_{n = 1}^{N} I_{n} (i, j) \sin (2 n π / N)}{\sum_{n = 1}^{N} I_{n} (i, j) \cos (2 n π / N)} \\ = atan \frac{\int \int_{x_{1 i}, y_{1 j}}^{x_{2 i}, y_{2 j}} \sum_{n = 1}^{N} I_{n} (x, y) \sin (2 n π / N) d x d y}{\int \int_{x_{1 i}, y_{1 j}}^{x_{2 i}, y_{2 j}} \sum_{n = 1}^{N} I_{n} (x, y) \cos (2 n π / N) d x d y} \end{matrix}

ϕ (i, j)

can be viewed as the argument of the vector

\vec{V} (i, j)

, Eq. (6).

Fig. 2 The schematic of CCD.

Download Full Size | PDF

\begin{matrix} \vec{V} (i, j) = [\int \int_{x_{1 i}, y_{1 j}}^{x_{2 i}, y_{2 j}} \sum_{n = 1}^{N} I_{n} (x, y) \cos (2 n π / N) d x d y] \vec{i} \\ + [\int \int_{x_{1 i}, y_{1 j}}^{x_{2 i}, y_{2 j}} \sum_{n = 1}^{N} I_{n} (x, y) \sin (2 n π / N) d x d y] \vec{j} \\ = \int \int_{x_{1 i}, y_{1 j}}^{x_{2 i}, y_{2 j}} [\sum_{n = 1}^{N} I_{n} (x, y) \cos (2 n π / N) \vec{i} + \sum_{n = 1}^{N} I_{n} (x, y) \sin (2 n π / N) \vec{j}] d x d y \\ = \int \int_{x_{1 i}, y_{1 j}}^{x_{2 i}, y_{2 j}} \vec{V} (x, y) d x d y \end{matrix}

It is clear that $\vec{V} (i, j)$ is the sum of all vectors in rectangle $(x_{1 i} : x_{2 i}, y_{1 j}, y_{2 j})$ . According to Eq. (1), the phase changes with the height. In practical measurement, the height is different at different point in the CCD cell. The wrapped phase $ϕ (i, j)$ , the argument of the vector sum, is not the true phase of any point in the cell. Therefore the errors are unavoidable because of the sampling, even under ideal condition.

We give an illustration, Fig. 3 . The height changes a little in the cell $(i_{1}, j_{1})$ . $ϕ (i_{1}, j_{1})$ , the argument of vector sum $\vec{V} (i_{1}, j_{1})$ , is nearly equal to $ϕ (x, y)$ , the argument of the vector $\vec{V} (x, y)$ [Fig. 4(a) ]. The error between $ϕ (i_{1}, j_{1})$ and $ϕ (x, y)$ is acceptable. This happens in most measurements because the CCD cell is small enough and the surface is smooth enough. But $ϕ (x, y)$ changes greatly in different part when the CCD cell is positioned at the abrupt discontinuities [Fig. 4(b)]. Equation (1) is inapplicable in this case because of the vector summing. This will cause the phase calculated from Eq. (5) become meaningless.

Fig. 3 the intensity distribution of CCD cells.

Download Full Size | PDF

Fig. 4 (a) The vector summing in cell (i ₁,j ₁) (b) the vector summing in cell (i ₂,j ₂).

Download Full Size | PDF

The area ratio of part A and B is $k_{A}$ and $\begin{matrix} k_{B} & (0 < k_{A}, k_{B} < 1, k_{A} + k_{B} = 1) \end{matrix}$ . We assume that the height of the object in part A and B changes a little separately. Intensity of part A and B are $I_{A n} (i_{2}, j_{2})$ and $I_{B n} (i_{2}, j_{2})$ . $ϕ_{A} (i_{2}, j_{2})$ is the argument of the vector ${\vec{V}}_{A} (i_{2}, j_{2})$ .

{\vec{V}}_{A} (i_{2}, j_{2}) = [\sum_{n = 1}^{N} I_{A n} (i_{2}, j_{2}) \cos (2 n π / N)] \vec{i} + [\sum_{n = 1}^{N} I_{A n} (i_{2}, j_{2}) \sin (2 n π / N)]

ϕ_{B} (i_{2}, j_{2})

is the argument of the vector

{\vec{V}}_{B} (i_{2}, j_{2})

.

{\vec{V}}_{B} (i_{2}, j_{2}) = [\sum_{n = 1}^{N} I_{B n} (i_{2}, j_{2}) \cos (2 n π / N)] \vec{i} + [\sum_{n = 1}^{N} I_{B n} (i_{2}, j_{2}) \sin (2 n π / N)]

Intensity of pixel $(i_{2}, j_{2})$ is

I_{n} (i_{2}, j_{2}) = k_{A} I_{A n} + k_{B} I_{B n}

$ϕ (i_{2}, j_{2})$ is the argument of the vector $\vec{V} (i_{2}, j_{2})$ , Eq. (10).

\begin{matrix} \vec{V} (i_{2}, j_{2}) = [\sum_{n = 1}^{N} I_{n} (i_{2}, j_{2}) \cos (2 n π / N)] \vec{i} + [\sum_{n = 1}^{N} I_{n} (i_{2}, j_{2}) \sin (2 n π / N)] \vec{j} \\ = k_{A} [\sum_{n = 1}^{N} I_{A n} (i_{2}, j_{2}) \cos (2 n π / N) \vec{i} + \sum_{n = 1}^{N} I_{A n} (i_{2}, j_{2}) \sin (2 n π / N) \vec{j}] \\ + k_{B} [\sum_{n = 1}^{N} I_{B n} (i_{2}, j_{2}) \cos (2 n π / N) \vec{i} + \sum_{n = 1}^{N} I_{B n} (i_{2}, j_{2}) \sin (2 n π / N) \vec{j}] \\ = k_{A} {\vec{V}}_{A} (i_{2}, j_{2}) + k_{B} {\vec{V}}_{B} (i_{2}, j_{2}) \end{matrix}

The expression above shows that vector $\vec{V} (i_{2}, j_{2})$ is the sum of vector $k {}_{A}{\vec{V}}_{A} (i_{2}, j_{2})$ and vector $k_{B} {\vec{V}}_{B} (i_{2}, j_{2})$ . We can encounter the following situations (Table 1 ).

Table 1. the results of vector summing

View Table

There are big errors in all cases. Especially in case 3 and 5, $ϕ (i_{2}, j_{2})$ is greater or less than $ϕ_{A} (i_{2}, j_{2})$ and $ϕ_{B} (i_{2}, j_{2})$ . There will be sharp spikes in the wrapped phases.

By the preceding analysis we know that CCD sampling will inevitably produce errors. The errors are acceptable at the smooth region on the object surface. Equation (1) is still applicable in such cases. The errors are dramatically big at the abrupt discontinuities. Equation (1) is no longer applicable in such cases. Therefore the abrupt discontinuities can be identified by observing whether the relationship between the phase and the equivalent wavelength complies Eq. (1) when the equivalent wavelength changes.

3. Invalid pixel identification

If Eq. (1) is applicable, the unwrapped phase is in inverse relation to the equivalent wavelength. The unwrapped phase linearly increases with the decrease of the equivalent wavelength in TPU. The gratings with difference equivalent wavelength (or spatial frequency) are projected onto the object one by one. If an exponentially growing sequence of grating spatial frequency is used, the least-squares fitting (LSF) is employed to get the estimator ω ^[12].

ω = \sum_{v = 0}^{\log_{2} s} 2^{v} φ_{u} (2^{v}) / \sum_{v = 0}^{\log_{2} s} 2^{2 v}

φ_{u} (2^{v})

is the unwrapped phase which is obtained by applying TPU. s is the number of gratings. LSF improves the accuracy of the measurement. ω is used to calculate the result of phase unwrapping. The standard deviation of LSF indicates the deviation between the data and the fitting curve. The standard deviation is small at the smooth region on the object surface. Equation (1) is inapplicable at the abrupt discontinuities because of the sampling. The standard deviation is dramatically big. If a suitable threshold of standard deviation is set, the abrupt discontinuities can be identified by comparing the standard deviation to the threshold. The pixels are positioned at the abrupt discontinuities if their standard deviations are bigger than the threshold. They are invalid. The phases of these pixels are deleted from the unwrapped phase map and recalculated by interpolating.

4. Simulation and experiment

4.1 Simulation demonstration

An object with both smooth surface and abrupt discontinuities is generated to demonstrate the theory in Matlab. The object is shown in Fig. 5 . $P_{0}$ takes 12.8cm, 6.4cm, 3.2cm, 1.6cm, and 0.8cm. System parameter θ is 27°. The deformed fringe patterns are shown in Fig. 6 . For a better demonstration of the effect of CCD sampling, both the shadow and noises are ignored. The height distribution obtained by temporal phase unwrapping is shown in Fig. 7 . The phases and the fitting curves at valid pixels (32, 36), (32, 38) and invalid pixel (32, 37) are shown in Fig. 8 . The true height values of pixel (32, 36) and (32, 38) are 0 cm and 0.5 cm separately. Pixel (32, 37) is positioned at the abrupt discontinuities. The reconstructed height at pixel (32, 37) is 0.59097cm. It is greater than its true value 0.5cm and that of pixel (32, 38). At the valid pixel (32, 36) and (32, 38), the measurements and fitting curves are consistent with the standard deviations 8.2516 × 10⁻¹⁵rad and 0rad separately. There is a big standard deviation 6.3842rad at the invalid pixel (32, 37).The mask which is shown in Fig. 9 is generated from standard deviation distribution using threshold 0.1rad. The white pixels in the mask are invalid pixels to be deleted and interpolated. The final result is shown in Fig. 10 .

Fig. 5 The object with abrupt discontinuities in simulation.

Download Full Size | PDF

Fig. 6 The deformed fringe patterns (P ₀ = 12.8cm, 6.4cm, 3.2cm, 1.6cm, 0.8cm, 0 θ = 27°).

Download Full Size | PDF

Fig. 7 The reconstructed object obtained by TPU.

Download Full Size | PDF

Fig. 8 Measurements and fitting curves.

Download Full Size | PDF

Fig. 9 The mask generated from standard deviations.

Download Full Size | PDF

Fig. 10 The reconstructed object after deleting and interpolation.

Download Full Size | PDF

4.2 Experiment

We measured a dental model with complex surface. The five-step phase-shifting is adopted in practical measurement. The deformed fringe patterns with equivalent wavelengths 32, 16, 8, 4, 2, 1 are shown in Fig. 11 . We excluded all pixels in the shadow and pixels with low modulation before using TPU. The reconstructed phase distribution of TPU is shown in Fig. 12 . It is apparent that big errors exist at the abrupt discontinuities. The standard deviations range from 0rad to about 30rad (Fig. 13 ). The standard deviations are much greater at the abrupt discontinuities than the smooth surface. We deleted all invalid points based on the mask shown in Fig. 14 and obtained the final phase distribution using two-dimensional linear interpolation (Fig. 15 ).

Fig. 11 The deformed fringe patterns.

Download Full Size | PDF

Fig. 12 The phase distribution obtained from TPU.

Download Full Size | PDF

Fig. 13 The grayscale of the standard deviations.

Download Full Size | PDF

Fig. 14 The mask generated from standard deviations.

Download Full Size | PDF

Fig. 15 The final phase distribution after deleting and interpolation.

Download Full Size | PDF

5. Conclusions

Many factors will affect the accuracy of PMP. The CCD sampling and the abrupt discontinuities on the object surface are two of them and not included in other algorithms. We rewrite the formula of phase retrieving and analyze the relationship between accuracy and CCD sampling by viewing the phase as the argument of a vector. The theoretical analysis shows that the CCD sampling introduces errors into the measurement and the phases obtained by N-step algorithm are not in inverse relation to the equivalent wavelength at the abrupt discontinuities. Therefore TPU is employed to identify the abrupt discontinuities. We presented the simulation demonstration and experiment. All results show the correctness of our theoretical analysis and validity of TPU in identification of phase errors caused by CCD sampling and the abrupt discontinuities on the object surface.

Acknowledgement

This project was supported by the National Natural Science Foundation of China (No 60838002).

References and links

1. V. Srinivasan, H. C. Liu, and M. Halioua, “Automated phase-measuring profilometry of 3-D diffuse objects,” Appl. Opt. 23(18), 3105–3108 (1984). [CrossRef] [PubMed]

2. J. M. Huntley, “Noise-immune phase unwrapping algorithm,” Appl. Opt. 28(16), 3268–3270 (1989). [CrossRef] [PubMed]

3. T. R. Judge and P. J. Bryanston-Cross, “A review of phase unwrapping techniques in fringe analysis,” Opt. Lasers Eng. 21(4), 199–239 (1994). [CrossRef]

4. H. Su, J. Li, and X. Su, “Phase algorithm without the influence of carrier frequency,” Opt. Eng. 36(6), 1799–1805 (1997). [CrossRef]

5. A. Asundi and Z. Wensen, “Fast phase-unwrapping algorithm based on a gray-scale mask and flood fill,” Appl. Opt. 37(23), 5416–5420 (1998). [CrossRef]

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

7. Y. Hao, Y. Zhao, and D. Li, “Multifrequency grating projection profilometry based on the nonlinear excess fraction method,” Appl. Opt. 38(19), 4106–4110 (1999). [CrossRef]

8. E. Hu, Y. He, and Y. Chen, “Study on a novel phase-recovering algorithm for partial intensity saturation in digital projection grating phase-shifting profilometry,” Opt. Int. J. Light Electron. Opt. (2008), doi:.

9. E. Hu, Y. He, and Y. Chen, “Study on a novel phase-recovering algorithm for partial intensity saturation in digital projection grating phase-shifting profilometry,” Optik - International Journal for Light and Electron Optics 121(1), 23–28 (2010). [CrossRef]

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

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

No.	The arguments of ${\vec{V}}_{A}$ and ${\vec{V}}_{B}$	The arguments of vector sum
1	$φ_{A} > 0, φ_{B} > 0, \| φ_{A} - φ_{B} \| < π$	$0 < φ_{A} < φ < φ_{B}$
2	$φ_{A} < 0, φ_{B} < 0, \| φ_{A} - φ_{B} \| < π$	$φ_{B} < φ < φ_{A} < 0$
3	$φ_{A} > 0, φ_{B} < 0, \| φ_{A} - φ_{B} \| > π$	$φ_{B} < 0 < φ_{A} < φ$
4	$φ_{A} > 0, φ_{B} < 0, \| φ_{A} - φ_{B} \| < π$	$φ_{B} < 0 < φ < φ_{A}$
5	$φ_{A} < 0, φ_{B} > 0, \| φ_{A} - φ_{B} \| > π$	$φ < 0 < φ_{A} < φ_{B}$
6	$φ_{A} < 0, φ_{B} > 0, \| φ_{A} - φ_{B} \| < π$	$φ_{A} < 0 < φ < φ_{B}$

Analysis and identification of phase error in phase measuring profilometry

Abstract

1. Introduction

2. Effect of CCD sampling on PMP

3. Invalid pixel identification

4. Simulation and experiment

4.1 Simulation demonstration

4.2 Experiment

5. Conclusions

Acknowledgement

References and links

Cited By

Figures (15)

Tables (1)

Equations (11)

Optics Express