Three-dimensional surface shape measurement based on grating projection is widely applied in various fields such as computer vision, physical simulation, automatic detection, biology, and medicine. It has many advantages, such as high speed, high accuracy, non-contact, automation, etc^[1-3]. Many scholars around the world have researched it and have achieved good results. For example, to reduce phase measurement errors caused by measuring glossy surfaces, Zhou et al. proposed a pixel-by-pixel combination multi-intensity matrix projection method. It has significantly reduced the number of projection operations and time consumption^[4]. Due to the influence of too many projected patterns on phase unwrapping, Yang et al. proposed a high-speed measurement method suitable for three-dimensional shapes which uses only three high-frequency internal shift phase modes (70 cycles), which improved measurement accuracy and speed, and obtained wrapping phase and fringe order^[5]. Using the gray code method, Lu et al. proposed a method based on staggered gray code light, which avoided step errors without projecting additional gray code modes^[6]. Due to overexposure in optical three-dimensional measurement, phase information cannot be obtained reliably. Feng et al. put forward a highly reflective surface measurement method based on pixel-by-pixel modulation, thereby significantly improving the measurement speed and accuracy^[7].

During the measurement process, due to the fact that the CCD (Charge Coupled Device) imaging system is not a general translation invariant system but a sampling imaging system with discrete characteristics, it will cause image distortion, resulting in spectra overlapping during the transformation process, and thus bringing errors to the measurement of the three-dimensional surface shape^[8-10]. Many scholars have proposed better measurement methods to reduce or eliminate these errors and improve accuracy^[11-13]. Due to the nonlinear effect of CCD on high-power laser wavefront detection, Du et al. proposed a method to reduce or eliminate spectra overlapping caused by the CCD nonlinear effect by increasing spatial carrier frequency^[11]. Due to the important role played by the nonlinearity of scientific-level CCD in experimental processing, Cheng et al. used two different schemes to experimentally test the nonlinear characteristics of CCD. They achieved good experimental results^[12]. Due to CCD’s nonlinear effects, spectra overlapping can occur when measuring complex optical three-dimensional surface shapes. Qiao et al. used the dual-frequency grating projection method to eliminate the nonlinear effects of CCD and improved measurement^[13].

The following sections study the influence of sampling on the measurement of three-dimensional surface shapes and provide detailed reasoning and analysis of its basic principles. They also effectively validate the basic principle analysis by conducting simulations and experiments that achieve improved results.

Fig. 1 shows the measurement system schematic diagram. $ {P_1}{P_2} $ is the projector’s optical axis, $ {L_0} $ is the distance between the optical center $ {I_2} $ of the CCD imaging system and the reference plane, $A$ and $C$ are the points located on the reference plane, and $D$ is the point on the object surface, $ h $ is the distance from $D$ to the reference surface.

By projecting the grating onto the surface of a three-dimensional object, the signal strength obtained by the CCD imaging system can be expressed as follows

$ \begin{split} g(x,y) = &r(x,y)\sum\limits_{n = - \infty }^\infty {{a_n}} \exp \{ j[2{\text{π}} n{f_0}x + n\phi (x,y)]\} = \\ & \sum\limits_{n = - \infty }^\infty {{q_n}} (x,y)\exp (j2{\text{π}} n{f_0}x)\quad, \end{split} $ (1)

where $r(x,y)$ is the non-uniform reflectivity of the object surface, $n$ is the Fourier series, ${a_n}$ is the n-order Fourier coefficient of $g(x,y)$, and $ {f_0} $ is the fundamental frequency of grating. And $\phi (x,y)$ is the phase of the object, where ${q_n}(x,y) = {a_n}r(x,y) \exp [jn\varphi (x,y)]$.

By performing the fast Fourier transform on Eq. (1) and using ${\text{π}} $ phase-shifting technology^[14] to eliminate the zero-order spectra components present in the frequency domain, the spectra expression containing the object height information can be obtained as follows

$ \begin{split} G({f_x},{f_y}) =& \int_{ - \infty }^\infty {\int_{ - \infty }^\infty {g(x,y)\exp [ - j2{\text{π}} ({f_x}x + {f_y}y)]{\mathrm{d}}x{\mathrm{d}}y} } = \\ & \sum\limits_{n = - \infty }^\infty {{Q_n}({f_x} - n{f_0},{f_y} - n{f_0})} +\\ &\sum\limits_{n = - \infty }^\infty {Q_n^*({f_x} - n{f_0},{f_y} - n{f_0})}\quad,\\[-8pt] \end{split} $ (2)

where $G({f_x},{f_{\text{y}}})$ and ${Q_n}({f_x},{f_y})$ are the spectra obtained by Fourier transform of $ g(x,y) $ and $ {q_n}(x,y) $, respectively, and $Q_n^*({f_x},{f_y})$ is the conjugate complex of ${Q_n}({f_x},{f_y})$.

Because CCD comprises an array of several small pixels arranged neatly and tightly with a certain geometric size, each CCD has approximately hundreds of thousands or even millions of pixels^[15]. Let the pixel shape be a rectangle, represented by the function $rect({x \mathord{\left/ {\vphantom {x {\Delta x}}} \right. } {\Delta x}},{y \mathord{\left/ {\vphantom {y {\Delta y}}} \right. } {\Delta y}})$, where $ \Delta x $ and $\Delta y$ present the pixel’s dimensions in the direction of the x-axis and y-axis, respectively. Therefore, the signal strength of each pixel can be represented as the convolution of $ g(x,y) $ and $rect({x \mathord{\left/ {\vphantom {x {\Delta x}}} \right. } {\Delta x}},{y \mathord{\left/ {\vphantom {y {\Delta y}}} \right. } {\Delta y}})$, its expression is as follows

$ g'(x,y) = g(x,y)*rect\left(\frac{x}{{\Delta x}},\frac{y}{{\Delta y}}\right) \quad.$ (3)

Using a comb function to sample Eq. (3), the discrete deformation fringe expression is obtained as follows

$ \begin{split} {g{''}}(x,y) =& {g{'}}(x,y)comb\left(\frac{x}{{\Delta {x_1}}},\frac{y}{{\Delta {y_1}}}\right) = \\ & [g(x,y)*rect\left(\frac{x}{{\Delta x}},\frac{y}{{\Delta y}}\right)]comb\left(\frac{x}{{\Delta {x_1}}},\frac{y}{{\Delta {y_1}}}\right), \end{split} $ (4)

where $ \Delta {x_1} $ and $ \Delta {y_1} $ present the sampling intervals of the fringes in the direction of the x-axis and y-axis, respectively.

The spectrum function obtained by performing Fourier transform on Eq. (4) is:

$ \begin{split} {G}''\left({f}_{x},{f}_{y}\right)=&G\left({f}_{x},{f}_{y}\right){\displaystyle \sum _{{n}_{x}=-\infty }^{+\infty }{\displaystyle \sum _{{n}_{y}=-\infty }^{+\infty }\delta \left({f}_{x}-\frac{{n}_{x}}{\Delta {x}_{1}},{f}_{y}-\frac{{n}_{y}}{\Delta {y}_{1}}\right)}}=\\ &\frac{1}{sl}{\displaystyle \sum _{{n}_{x}=-\infty }^{+\infty }{\displaystyle \sum _{{n}_{y}=-\infty }^{+\infty }\mathrm{sin}c\left(\frac{{n}_{x}{\text{π}} }{s}\right)\mathrm{sin}c\left(\frac{{n}_{y}{\text{π}} }{l}\right)Q\left({f}_{x}-{f}_{0}-\frac{{n}_{x}}{\Delta {x}_{1}},{f}_{y}-{f}_{0}-\frac{{n}_{y}}{\Delta {y}_{1}}\right)}}\quad,\end{split} $ (5)

where $ {n_x} $ and ${n_y}$ present pixel points on the x-axis and y-axis, respectively, $s = {{\Delta {x_1}} \mathord{\left/ {\vphantom {{\Delta {x_1}} {(\Delta x + \Delta {x{'}})}}} \right. } {(\Delta x + \Delta {x{'}})}}$ and $l = {{\Delta {y_1}} \mathord{\left/ {\vphantom {{\Delta {y_1}} {(\Delta y + \Delta y')}}} \right. } {(\Delta y + \Delta y')}}$ present the number of sampling points for each fringe in the corresponding direction, respectively, and $ \Delta {x{'}} $ and $ \Delta {y{'}} $ present the pixel spacing in the corresponding direction, respectively.

Thus, the spectra of the sampling function are the infinite repetition for the spectra of the primitive continuous function in the frequency domain. This is commonly known as "spectra island"^[15-19]. As a result, in addition to $ {f_0} $, the higher-order spectral components, such as second and third order, are also generated.

Due to the useful information containing changes in object height within $ {f_0} $ of the spectra, a suitable low-pass filter must be designed to gain the $ {f_0} $ and remove the high-order spectra component.

The low-pass filter is a modulation system of a point spread function. Its filtering process is the convolution process of the spectrum function and the point spread function $ {G_1}({f_x},{f_y}) $. $ {G_1}({f_x},{f_y}) $ in the frequency domain is presented by a Gaussian filter as follows

$ {G_1}({f_x},{f_y}) = \frac{1}{{2{\text{π}} }}\exp \left( - \frac{{{f_x}^2 + {f_y}^2}}{{2{\delta ^2}}}\right) \quad,$ (6)

where $ \delta $ presents the standard deviation of the filter related to the degree of defocus.

The spectra signal obtained through system defocusing is

$ {G{''}}{'}({f_x},{f_y}) = {G{''}}({f_x},{f_y}) \otimes {G_1}({f_x},{f_y})\quad. $ (7)

By filtering, the higher-order harmonics can be well separated from $ {f_0} $, so that only one of the $({1 / {sl}})\sin c({{\text{π}} /s})\sin c({{\text{π}}/ l})Q({f_x} - {f_0} - {1 / {\Delta {x_1}}},{f_y} - {f_0} - {1 /{\Delta {y_1}}})$ is retained after filtering out the higher-order harmonic components.

Then, the inverse Fourier transform is enforced on Eq. (7) to reconstruct the signal strength taken by the CCD imaging system. The measurement system outputs the n-th sine fringe image, which is obtained by the CCD system, as shown below

$ \begin{split} g^{\land}(x,y) =& {F^{ - 1}}[{G{''}}{'}({f_x},{f_y})] = \\ &\sum\limits_{k = {\text{ - }}\infty }^\infty {A_k^{\land}\cos \{ k[2{\text{π}} {f_0}x + \phi (x,y) + {\delta _n}]\} }, \end{split} $ (8)

where $ {F^{ - 1}}[ \cdot ] $ presents inverse Fourier transform, $ A_k^{\land} $ presents the Fourier coefficient of $ g_n^{\land}(x,y) $, and $ {\delta _n} $ presents the phase-shift amount, $ {\delta _n}{\text{ = }}{{2n{\text{π}} } /{{n_1}}} $, $ n = 1,2, \cdots , {n_1} $.

When the n-step phase-shift method is used, we can gain the phase as follows

$ {\phi ^ \wedge }(x,y) = {\mathrm{arctan}}\left[\frac{{\displaystyle\sum_{n = 1}^N {g^ \wedge _n(x,y)} \sin ({\delta _n})}}{{\displaystyle\sum_{n = 1}^N {g^ \wedge _n(x,y)} \cos ({\delta _n})}}\right]\quad, $ (9)

where $ {\phi ^ \wedge }(x,y) $ is the wrapped phase.

Under the conditions of the telecentric projection optical path, considering ${L_o} \gg h(x,y)$ in the real measurement conditions, $h(x,y)$ and $ {\phi ^ \wedge }(x,y) $ will satisfy the relations

$ h(x,y) = - \frac{{L{\phi ^ \wedge}(x,y)}}{{2{\text{π}} {f_0}d}}\quad.$ (10)

It has been discussed in the literature [15] that in order to ensure the separation of the $ {f_0} $ from other periodic spectra components and the separation of spectra components during the same period, the sampling condition $ m > 4 $ (where $m = {{\Delta f} \mathord{\left/ {\vphantom {{\Delta f} {{f_0}}}} \right. } {{f_0}}}$, $\Delta f$ presents sampling frequency) must be satisfied. This can avoid overlapping between the $ {f_0} $ and the higher-order spectra components, so as to accurately reconstruct the object shape measured. Otherwise, reconstruction is difficult.

Lastly, according to $s = {{\Delta {x_1}} \mathord{\left/ {\vphantom {{\Delta {x_1}} {{\text{(}}\Delta x + \Delta {x{'}}{\text{)}}}}} \right. } {{\text{(}}\Delta x + \Delta {x{'}}{\text{)}}}}$ and $l = {{\Delta {y_1}} \mathord{\left/ {\vphantom {{\Delta {y_1}} {(\Delta y + \Delta {y{'}})}}} \right. } {(\Delta y + \Delta {y{'}})}}$, combining the relationship $\Delta {f_x} = {1 \mathord{\left/ {\vphantom {1 {\Delta {x_1}}}} \right. } {\Delta {x_1}}}$ as well as $\Delta {f_y} = {1 \mathord{\left/ {\vphantom {1 {\Delta y}}} \right. } {\Delta y}}$ between sampling frequency and sampling interval, it can be obtained that

$ \left\{ \begin{gathered} {m_x} = \frac{1}{{\Delta {x_1}{f_0}}} = \frac{1}{{s(\Delta x + \Delta {x{'}}){f_0}}} \\ {m_y} = \frac{1}{{\Delta {y_1}{f_0}}} = \frac{1}{{l(\Delta y + \Delta {y{'}}){f_0}}} \\ \end{gathered} \right. \quad,$ (11)

where ${m_x}$ and ${m_y}$ present the sampling frequency ratio in the direction of x-axis and y-axis to ${f_0}$, respectively.

It can be seen that the method of reducing the sampling interval, i.e., the number of sampling points per fringe, to increase $ m $, can be used to increase the accuracy of object surface shape measurements.

Simulation and experiment were executed to validate the basic principle analysis.

We performed computer simulation verification on the analysis of basic principles. Assuming that the geometric parameter of the measurement system is $ {{{L_0}} \mathord{\left/ {\vphantom {{{L_0}} d}} \right. } d} = 4 $. The simulated object surface shape is shown in Fig. 2 (color online), with a size of 512×512 pixels.

We projected a digital projector onto a simulated object and used a CCD camera system to obtain deformation fringes. If 40 fringes are taken, then ${f_0} = {{40} \mathord{\left/ {\vphantom {{40} {512}}} \right. } {512}}$ fringe/pixel. Using MATLAB to process the fringes, the sampling intervals of both the x-axis and y-axis directions were 8 pixels, thus, $m = 1.600\;0$ can be obtained from Eq. (11). Reducing the sampling interval of the fringes to 0.5 times the original sampling interval, i.e., the sampling interval was 4 pixels. From $s = {{\Delta {x_1}} \mathord{\left/ {\vphantom {{\Delta {x_1}} {{\text{(}}\Delta x + \Delta {x{'}}{\text{)}}}}} \right. } {{\text{(}}\Delta x + \Delta {x{'}}{\text{)}}}}$ and $l = {{\Delta {y_1}} \mathord{\left/ {\vphantom {{\Delta {y_1}} {(\Delta y + \Delta {y{'}})}}} \right. } {(\Delta y + \Delta {y{'}})}}$, it can be seen that the number of sampling points for each fringe was 0.5 times that of the original, and from Eq. (11), $m = 3.200\;0$ can be obtained. It can be seen that neither of the two situations satisfied the sampling condition $ m > 4 $. The obtained spectra diagrams are shown in Fig. 3 (a) and 3(b) (color online), respectively.

Then, we made the sampling points of each fringe ${1 \mathord{\left/ {\vphantom {1 4}} \right. } 4}$ and ${1 \mathord{\left/ {\vphantom {1 8}} \right. } 8}$ times the original, i.e., the sampling interval was 2 pixels and 1 pixel, respectively. $m = 6.400\;0$ and $m = 12.800\;0$ can be obtained, respectively, indicating that the sampling condition $ m > 4 $ was satisfied. The obtained spectra diagrams are shown in Fig. 3 (c) and 3(d) (color online), respectively.

In Figs. 3 (a) and 3(b), the components of ${f_0}$ in the spectra diagrams overlap with the higher-step frequency components because the sampling condition $ m > 4 $ was not satisfied. But in Figs. 3 (c) and 3(d), the corresponding frequency components are separated because the sampling condition was satisfied. In the four sub-figures, the smaller the number of sampling points, the larger the $m$, and the better the separation effect.

Fig. 4 shows the surface shape errors between the reconstructed and simulated objects in the above four situations.

The maximum absolute error value (MAEV) and average absolute error value (AAEV) of each sub-image in Fig. 4 are shown in Tab. 1. The MAEV obtained from the last three sampling intervals are 88.80%, 38.38%, and 31.50%, and the AAEV are 71.84%, 43.27%, and 32.26%, respectively, of the first sampling interval.

As shown by Fig. 4 and Tab. 1, when sampling condition $ m > 4 $ is not satisfied, it is difficult to reconstruct the object’s surface shape and has relatively large errors. On the contrary, the error is relatively small. The larger $m$, the smaller the error, but improving the resolution of the CCD imaging system is necessary.

To further validate the impact of CCD imaging system sampling on measuring three-dimensional surface shapes, an actual measurement experiment of a hemispherical object was carried out. As shown in Fig. 5, the simple experimental device system uses a digital projector and a low-distortion CCD camera.

The actual experiment used the same method as the computer simulation to obtain the deformation fringes of the experimental object. 600×480 pixels and 40 fringes were taken, the fundamental frequencies were ${f_x} = {{40} \mathord{\left/ {\vphantom {{40} {600}}} \right. } {600}}$ fringe/pixel and ${f_y} = {{40} \mathord{\left/ {\vphantom {{40} {480}}} \right. } {480}}$ fringe/pixel, respectively. Using MATLAB to process the fringes, the sampling intervals of two directions were 8 pixels. From Eq. (11), ${m_x} = 1.875\;0$ and ${m_y} = 1.500\;0$ can be obtained.

Using the same method as computer simulation, the fringes’ sampling interval was reduced to 0.5 times that of the original fringes, changing it to 4 pixels. Similarly, ${m_x} = 3.750\;0$ and ${m_y} = 3.000\;0$ can be obtained from Eq. (11).

Neither of the above situations satisfies the sampling condition $ m > 4 $. The reconstructed object surface shapes are shown in Fig. 6 (a) and 6(b) (color online).

When the sampling intervals for each fringe in the directions of the x-axis and y-axis are set to 2 pixels and 1 pixel, respectively, then ${m_x} = 5.625\;0$, ${m_y} = 4.500\;0$, and ${m_x} = 7.500\;0$, ${m_y} = 6.000\;0$ can be obtained, respectively. Both situations satisfy the sampling condition $ m > 4 $. The reconstructed object surface shapes are shown in Fig. 6 (c) and 6(d) (color online).

Comparing the experimental reconstruction results of the four sub-images in Fig. 6 yields the same conclusions as the computer simulation results mentioned above.

Due to the impact of sampling on the accuracy of three-dimensional surface shape measurements, the basic principles of spectra overlapping, spectra separation, and measurement accuracy caused by sampling were analyzed and discussed.

The comb function was used to sample the signal intensity of each pixel in the CCD, and discrete deformation fringes were obtained. After the Fourier transform of the deformation fringes, the "spectra island" is generated in the frequency domain. When the method of reducing the CCD sampling interval is used, i.e., the sampling frequency is increased to satisfy the sampling condition $ m > 4 $, then the adjacent " spectra islands" will not overlap with the ${f_0}$; thus, reconstruction of the three-dimensional shape can be improved.

Reducing the sampling interval of pixels causes the ratio $m$ to increase. When the sampling conditions are satisfied, the larger the $m$, the better the recovery of the object’s three-dimensional shape. However, further improvement in the resolution of the CCD imaging system is necessary.

In the computer simulation and practical experiment, the smaller the sampling interval for each fringe, i.e., the fewer sampling points, the better the reconstruction of the object’s surface shape. When sampling condition $ m > 4 $ is satisfied, it effectively reconstructs the object’s three-dimensional shape.