## Abstract

In this paper, we propose a phase difference minimization algorithm to measure the specular surface shape in a displacement-free stereoscopic phase measuring deflectometry (PMD) system. The presented system is capable of solving the height-normal ambiguity appearing in a PMD system without moving any system component. Both the surface normal and the absolute height are simultaneously obtained by implementing phase difference minimization between the phase distributions in the LCD screen and the camera image plane. In particular, phase difference minimization is performed by using a second order polynomial fitting iteration method. Bi-cubic sub-pixel interpolation combined with 2D Fourier integration is used to reconstruct the surface. Finally, the performance of the proposed stereoscopic PMD system is verified by measuring the surface shapes of different mirrors and performing repeatability tests.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Specular surface shape measuring techniques have been excessively investigated in recent years. Among these techniques, interferometry and phase measuring deflectometry demonstrate great advantages to obtain three-dimensional surface shapes. An interferometer guarantees high precision metrology of regular shape specular surfaces (i.e., flat, spherical or aspherical mirrors) using an extremely high precision reference mirror [1], but it finds difficulty to measure the shapes of complicated specular surfaces such as large curvature radius mirrors or free-form mirrors. Hence, phase measuring deflectometry (PMD) technique [2], which is able to obtain both the surface height and normal of a specular surface, is proposed to complement the existing measurement methods. Compared to other specular surface measurement systems, PMD demonstrates considerable advantages as it guarantees both a large measuring area and an easy implementation by maintaining high measuring precision. What is more, it is even able to measure the surface shapes of those specular objects which contain isolated regions [3].

PMD measurement principle relies on the very simple “Law of reflection” from which the surface shape is determined by performing an inverse ray tracing from the camera image plane to the LCD screen plane. Unfortunately, the height-normal ambiguity, demonstrated as the surface normal and height cannot be determined simultaneously, is unavoidable in PMD measurement [2]. Thus, to solve this ambiguity, different techniques have been proposed. One classical solution to eliminate the height-normal ambiguity is to either move the screen or the camera to different positions to obtain a solitary reflection relation between the LCD panel and the camera image plane [4,5]. However, such physical displacement deteriorates the measurement precision if we consider component displacement errors. To avoid system component displacements, Burge *et al.* has proposed SCOTS system with the camera and the LCD screen fixed [6–8]. Even though SCOTS guarantees high accuracy, its application could be limited because SCOTS is more appropriate to measure a regular shape specular surface (i.e., spherical or aspherical mirrors), but it may not be suitable to measure a free-form surface. On the other hand, another PMD system has been proposed by applying monoscopic fringe reflectometry technique [9–11]. In this case, the height-normal ambiguity is regularized by introducing a prior surface height assumption. Afterwards, Huang *et al.* further improved monoscopic PMD by using a well-established mathematical model to assume the tested surface shape [12]. Then, by properly selecting a model to define the surface shape and later implementing the coefficient optimization with the captured fringe patterns, the surface can be reconstructed. Although the above-mentioned monoscopic systems [9–12] fulfill the specular surface shape measurement, the definition of the hypothetical reference specimen directly influences the measurement result, and thus high accuracy may not be guaranteed if the reference is not properly selected. Considering that limitation of monoscopic PMD, other PMD systems employing more screens or cameras are proposed. In particular, Zhang *et al.* [3,13] proposed an interesting PMD scheme by using two LCD screens with a beam-splitter. In their system, a beam-splitter images the fringe pattern on one screen to a certain plane parallel to the second LCD screen so these two fringe patterns (i.e., the real pattern on the second screen and the imaged pattern of the first screen) share a spatial distance and therefore, a virtual screen displacement is guaranteed by avoiding any physical movement. Similar to Zhang system, Li *et al.* as well suggested a bi-plane PMD scheme, where the bi-plane is obtained by directly covering a transparent panel above the LCD screen [14]. The above-mentioned systems guarantee an easy implementation and they provide the desirable robustness, but the final measurement precision is dramatically influenced by the surface homogeneity of the screen-size beam-splitter or transparent panel. Apart from adding one screen, an easier PMD scheme is proposed by Knauer *et al.* [2] where two cameras are introduced to capture the reflected fringe patterns. In their method, one point ** A** is selected in the first camera image plane and then it is linked to the optical center of the first camera

**to forge an inverse ray tracing. Hence, a reflection vector $\overrightarrow {{\boldsymbol {AO}}} $ is obtained and therefore, any potential surface height and normal corresponding to point**

*O***can be calculated along this vector. Subsequently, by assuming one point along the reflection vector $\overrightarrow {{\boldsymbol {AO}}} $ as the surface point, we connect it to the second camera optical center and it intersects to the second camera image plane at**

*A***, and we obtain the phase value of**

*B***. Next, we find a point in the LCD plane which has the same phase value of**

*B***, and we forge a secondary reflection relation between the second camera and the LCD by using the same phase value pair. At this time, the normal to the chosen point with respect to the second camera is also determined. Finally, the real surface point is obtained by searching along the whole reflection vector $\overrightarrow {{\boldsymbol {AO}}} $ until the two normal corresponding to both the first camera and the second camera coincide. The Knauer scheme demonstrates great advantages for its displacement-free characteristic, and it directly calculates the surface normal to the test object. Therefore, by accurately estimating an initial surface point which is close to the real surface, their method can achieve rapid surface point calculation [2]. Nevertheless, as it is usually not easy to directly assume the real surface point, a secondary surface point assumption is always required. The secondary surface point assumption is challenging because, without a clear indication, we may select this secondary point even farther to the initially assumed point, and this will increase the calculation time.**

*B*Considering the shortcomings of the above-discussed PMD systems, we proposed a new phase difference minimization algorithm in a stereoscopic displacement-free PMD system to perform specular surface shape measurement. Note that even though stereoscopic scheme has been broadly used in profilometry [15,16], it has not been widely implemented in deflectometry. In our system, the height-normal ambiguity is eliminated by introducing a second camera and thus, both the surface height and normal can be obtained simultaneously. Specifically, compared to Knauer system [2], our system uses a phase difference minimization method with a second order polynomial fitting iteration algorithm. As iteration process allows the calculated surface point to rapidly approach the real surface, our method allows a flexible initial surface point assumption, and it also promises short calculation time. On the other hand, our proposed stereoscopic PMD scheme can guarantee the system robustness as no physical displacement is required. To be more specific, in order to achieve height-normal ambiguity elimination and calculate the normal to the surface, we implement our phase difference minimization method by using the phase patterns in the second camera and the LCD. Finally, the surface is reconstructed using 2D Fourier integration [17].

The outline of this paper is as follows. In Section 2, we first briefly introduce height-normal ambiguity in a phase measuring deflectometry (PMD) system. Subsequently, we present our stereoscopic displacement-free PMD system and phase difference minimization principle for height-normal ambiguity elimination. In particular, we comprehensively discuss the implementation of second order polynomial fitting iteration for phase difference minimization. Moreover, we also present specular surface reconstruction in Section 2.3. The simulation of a flat specular surface measurement using our stereoscopic PMD system is demonstrated in Section 3. Specifically, we perform the simulation with both calibration-error-free and calibration-error-contained stereoscopic PMD systems. In this case, we have comprehensively studied the influence of calibration errors to surface measurement results. In Section 4, we calibrate a stereoscopic PMD system, and we measure the surface shapes of both a flat and a spherical mirror. Moreover, to further verify the system robustness, repeatability tests are performed. Finally, we provide the main conclusions in Section 5.

## 2. Stereoscopic displacement-free phase measuring deflectometry

#### 2.1 Height-normal ambiguity

Different from profilometry which measures diffuse objects via diffuse reflection with a video projector [18], the principle of deflectometry for specular surface shape measurement relies on the very basic “Law of reflection”. A digital screen (i.e., LCD) is introduced to display fringe patterns and then it is used to illuminate the surface under test. On the other side, we place a camera at a certain angle where it can capture the distorted fringe pattern image reflected by the specular surface. Here we want to emphasize that the distortion of fringe pattern between the LCD and the camera is caused by irregularity of the reflecting surface shape. The basic PMD scheme is shown in Fig. 1(a). The pattern on the LCD is a gray-level based sinusoidal fringe pattern with its intensity being modulated by a phase distribution. Then, we select an LCD pixel of phase value as *φ* [see point ** A’** in Fig. 1(b)], and this LCD pixel is reflected into one pixel in the camera image plane [point

**in Fig. 1(b)]. Moreover, as surface reflection only changes the intensity but not the phase value, these two pixels have the same phase value. Next, by regarding the camera as a pinhole camera [19], we can implement inverse ray tracing from the pinhole camera to the LCD by using this pixel pair, and the surface normal as well as its height are determined. Such inverse ray tracing process can be presented by the red vectors shown in Fig. 1(b), and we finally determine one surface point at**

*A***with its normal as $\overrightarrow {{{\boldsymbol n}_1}} $. In Fig. 1(b),**

*P***is the principle point of the pinhole camera (this is obtained by camera calibration [19]), $\overrightarrow {{\boldsymbol V}\; } $ is the inverse ray tracing vector from**

*O*_{c}**to**

*A***passing**

*P***.**

*O*_{c}The above-mentioned method seems promising to obtain the surface height and normal once we perform inverse ray tracing to all pixels in the pinhole camera. Unfortunately, the height-normal ambiguity, demonstrated as the surface normal and height cannot be uniquely determined, is unavoidable within this PMD scheme. Specifically, this can be seen in Fig. 1(b) by prolonging vector $\vec{{\boldsymbol V}}$ until it meets another point ** Q** on the second surface. By now, we can establish a secondary reflection relation as the blue rays in Fig. 1(b) and it gives a different height and normal (i.e., $\overrightarrow {{{\boldsymbol n}_2}} $). Therefore, the surface point determination in the basic PMD system is not available as we cannot decide which point (i.e.,

**or**

*P***) represents the real surface.**

*Q*#### 2.2 Stereoscopic phase measuring deflectometry

As height-normal ambiguity restricts us from finding the specular surface point, we use stereoscopic phase measuring deflectometry. Specifically, a stereoscopic PMD system consists of a flat LCD, two cameras and the test surface. The LCD displays sinusoidal fringe patterns (SFPs) while both cameras, located at proper positions, capture the reflected SFPs from the test surface. The intensity distribution of an SFP on the LCD is determined by a certain phase pattern *Φ _{L}*(

*x,y*). Subsequently, as the test surface distorts the SFP patterns, the phase pattern in any camera is represented as

*Φ*(

_{C}*i,j*). To establish the correspondence between the LCD and the cameras via phase values, we perform phase shifting [20] to the SFPs in the cameras, and we get the wrapped phase pattern. Next, we use temporal phase unwrapping [21,22] to further obtain the unwrapped phase distribution

*Φ*(

_{C}*i,j*). Note that as we focus on height-normal ambiguity elimination and surface measurement with stereoscopic PMD in this paper, we directly use the multi-frequency temporal phase unwrapping method described in Ref. [22].

Once we obtain the unwrapped phase distributions in both the cameras and the LCD, we start inverse ray tracing for surface point determination. First, we select any given pixel ** A_{1}** in the first camera image plane which has a phase value of

*Φ*, and we build a basic vector $\overrightarrow {{{\boldsymbol V}_1}} $ by connecting

_{A}**with the camera principle point**

*A*_{1}**. Here, the cameras are regarded as pinhole cameras, and this vector is presented in Fig. 2 as the red solid line. Next, we find on the LCD screen a certain pixel**

*O*_{C1}**which has a same phase value of**

*L*_{A}*Φ*. Now, because phase value is not changed via surface reflection, we can tell that the ray exiting from

_{A}**is reflected by the tested surface and then enters**

*L*_{A}**through**

*O*_{C1}**. Thus, the real surface point corresponding to**

*A*_{1}**must be located at one point along vector $\overrightarrow {{{\boldsymbol V}_1}} $. Unfortunately, height-normal ambiguity restricts us from finding this real surface point. Nevertheless, we can assume any point**

*A*_{1}**along vector $\overrightarrow {{{\boldsymbol V}_1}} $ as the surface point, and we connect**

*P***to**

*P***to forge vector $\vec{{\boldsymbol R}}$. At this moment, the surface normal to point**

*L*_{A}**can be calculated by Law of reflection because we already know both $\overrightarrow {{{\boldsymbol V}_1}} $ and $\vec{{\boldsymbol R}}$. For instance, the normal to three potential surface points (i.e.,**

*P***,**

*P*_{1}**and**

*P*_{2}**in Fig. 2) along vector $\overrightarrow {{{\boldsymbol V}_1}} $ are acquired by forging their inverse ray tracing. Here, we assume that**

*P*_{3}**is the real surface point, whereas**

*P*_{1}**and**

*P*_{2}**represent fake surface points.**

*P*_{3}For the next step, we introduce the second pinhole camera to find the real surface point. Here, let us again assume a point ** P** along vector $\overrightarrow {{{\boldsymbol V}_1}} $ is the surface point and then we intersect

**with the principle point of the second pinhole camera (**

*P***in Fig. 2) to forge another vector $\overrightarrow {{{\boldsymbol V}_2}} $. The intersection of $\overrightarrow {{{\boldsymbol V}_2}} {\; }$ to the second camera image plane gives us point**

*O*_{C2}**, and its phase value is**

*B**Φ*. Then, as the normal to the hypothetical surface point

_{B}**is calculated as $\vec{{\boldsymbol n}}$ by the first camera, we use $\vec{{\boldsymbol n}}$ to determine the inverse reflection point**

*P***on the LCD plane corresponding to**

*L*_{B}**, and the phase value of point**

*B***is**

*L*_{B}*Φ*.

_{LB}Now, if the assumed point ** P** we picked is the real surface point

**(**

*P*_{1}**coincides with**

*P***), then we obtain**

*P*_{1}**on the second camera image plane and**

*B*_{1}**on the LCD plane (see the green solid lines in Fig. 2). Moreover, the phase values**

*L*_{B1}*Φ*and

_{B}*Φ*must be equal, because these two points are indeed obtained from the real reflection. However, if the assumed point

_{LB}**we picked is the fake surface point**

*P***, then the corresponding points**

*P*_{2}**and**

*B*_{2}**(see yellow solid lines in Fig. 2) do not have a same phase value, so the phase difference is not zero. The reason to such phase difference is because in the real implementation,**

*L*_{B2}**is reflected by another surface point rather than**

*L*_{B2}**, and this causes the real reflected point to deviate from**

*P*_{2}**. Finally, we also analyze the last case as the assumed point is selected as**

*B*_{2}**(blue solid lines in Fig. 2). By comparing**

*P*_{3}**to**

*P*_{3}**, point**

*P*_{2}**is located farther from the real surface point**

*P*_{3}**along vector $\overrightarrow {{{\boldsymbol V}_1}} $. Thus, the phase difference between**

*P*_{1}**and**

*P*_{3}**should be greater than the**

*L*_{B3}**case.**

*P*_{2}From the above analysis, surface point determination can be easily performed by minimizing the absolute phase difference between the LCD and the second camera with respect to point ** P** along vector $\overrightarrow {{{\boldsymbol V}_1}} $. Moreover, we calculate the absolute phase differences of all assumed surface points along vector $\overrightarrow {{{\boldsymbol V}_1}} $, and we consider these absolute phase differences as a function. We surmise it is a smooth function with only one minimum, because we are assuming a smooth surface. Thus, by finding the minimum in this function, we obtain the real surface point. Here, we calculate the square of the phase difference as δ in Eq. (1),

*Φ*and

_{B}*Φ*are the corresponding unwrapped phase values on the second camera image plane and the LCD along vector $\overrightarrow {{{\boldsymbol V}_1}} $. In order to utilize Eq. (1) for phase difference minimization, we define another tiny value ε. Now, if δ of any assumed surface point is smaller than ε, we take this point as the real surface point. To accurately establish phase correspondence between the second camera and the LCD, we let the square root of ε be smaller than the phase difference of any two adjacent points in the unwrapped phase map. A smaller ε allows more iterations so the measurement precision is enhanced, but it also decreases the measurement speed. Ultimately, the value of ε is determined by the noise in the measurement. To compromise the measurement precision and speed, we selected ε of 0.01 in our experiments (see Section 4).

_{LB}We propose a second order polynomial fitting iteration method for finding the surface point by phase difference minimization. A flow chart diagram is shown in Fig. 3 for better visualization of this method. For the first step (see Step 1 in Fig. 3), we select several assumed surface points along vector $\overrightarrow {{{\boldsymbol V}_1}} $, and we calculate their phase difference squares [i.e. Equation (1)] by performing inverse ray tracing between the LCD and the second camera. Specifically, to implement a robust polynomial fitting, these assumed surface points should be selected close to the real measured surface. Next, in Step 2, a second order polynomial function is used to fit these phase difference values, and the minimum in this function is found where its phase difference square value is δ. Now, if δ is smaller than ε, we direct output this minimum point as the real surface point. Otherwise, if δ is greater than ε, we trigger a subsequent polynomial fitting iteration until δ becomes smaller than ε (see Step 4 to 6 in Fig. 3). Make notice that the assumed surface points vary within each iteration loop, and they are gradually approaching the real surface point. By properly assuming potential surface points, the runtime to terminate the iteration and output a measuring point in our experiments (see Section 4) is nearly 0.2 second using Matlab. As we have not made an optimization of the calculation time, it can be further improved by using parallel computing, as well as using the surface point of a previous pixel.

#### 2.3 Surface reconstruction

To obtain the surface normal with high accuracy for surface reconstruction, we further introduce sub-pixel interpolation. A pixel *P _{CCD}* in the first camera image plane with a phase value of

*Φ*is initially selected and we find its corresponding pixel in the LCD as

*P*(see these two pixels labeled red in Fig. 4). Here, the coordinates of both

_{LCD}*P*and

_{CCD}*P*need to be determined in order to accurately perform the inverse ray tracing. Thus, we use the geometric center of

_{LCD}*P*to represent its coordinate (see point

_{CCD}**in Fig. 4), but we cannot use the geometric center of**

*A**P*(i.e., point

_{LCD}**in Fig. 4) as its corresponding coordinate. This is because the LCD pixel finite size introduces a quantization phase error [23], so the phase value at**

*A*’**may not be equal to the phase value of**

*A*’**. Hence, we use bi-cubic interpolation to find another point**

*A***in pixel**

*A*’’*P*which has the phase value of

_{LCD}*Φ*for correct inverse ray tracing. To be more specific, we describe the sub-pixel coordinate of point

**by Eqs. (2) and (3),**

*A*’’*x*(

*Φ*,

_{hLCD}*Φ*) and

_{vLCD}*y*(

*Φ*,

_{hLCD}*Φ*) represent the sub-pixel coordinate of point

_{vLCD}**,**

*A*’’*Φ*and

_{hLCD}*Φ*are the horizontal and vertical phase values of any pixel,

_{vLCD}*a*and

_{xij}*a*are bi-cubic interpolation coefficients. Here, we have both horizontal and vertical phase values because the sinusoidal fringe patterns used in our system are presented in these two directions. In order to perform interpolation, we collect the phase values of both the studied pixel and its surrounding pixels, and we use the method presented in Ref. [24] to calculate the bi-cubic interpolation coefficients. Once the coefficients are obtained, we can use Eqs. (2) and (3) to easily calculate the coordinate of

_{yij}**in any LCD pixel.**

*A*’’Note that the above discussed bi-cubic sub-pixel interpolation is feasible to establish the correspondence between the first camera and the LCD. On the other hand, we also need to perform a secondary sub-pixel interpolation between the second camera and the LCD. Same bi-cubic interpolation principle is directly used for the secondary sub-pixel interpolation as it also depends on inverse ray tracing.

At this moment, by using phase minimization principle described in Section 2.2 and bi-cubic sub-pixel interpolation, the normal and coordinate of any surface points are determined. However, as the cameras and the LCD are slightly separated to maintain a large measurement area, the coordinate in *z* direction encounters larger ray tracing error than the (*x*,*y*) coordinates. Therefore, we introduce 2D Fourier integration method to reconstruct the surface. To do so, we entitle the surface shape as *f*(*x*,*y*), and the surface derivatives of any point are expressed as *f ^{x}*(

*x*,

*y*) and

*f*(

^{y}*x*,

*y*) in

*x*and

*y*directions, which can be easily calculated from the surface normal. In the real implementation, as we obtain the normal to the whole surface,

*f*(

^{x}*x*,

*y*) and

*f*(

^{y}*x*,

*y*) are in fact two matrices. Nevertheless, experimental measurements may give us two inhomogeneous derivative distributions, as the distance between any two adjacent pixels are not equal. So, we perform interpolation to sample them into homogeneously distributed matrices which are required by 2D Fourier integration. Once the derivative matrices are homogeneous, we transform the derivatives to the Fourier domain, and we obtain Eqs. (4) and (5),

We also perform Fourier transform to the second order derivatives, and we get Eqs. (6) and (7),

Now, we introduce Laplacian operator to the surface shape *f*(*x*,*y*), and we again perform Fourier transform to obtain Eq. (8),

Finally, by merging Eq. (4) to Eq. (7) into Eq. (8), we acquire the surface shape in the Fourier domain as:

*u*and

*ν*are the surface point coordinates in the Fourier domain,

*F*(

^{u}*u*,

*ν*) and

*F*(

^{v}*u*,

*ν*) are the surface derivatives in the Fourier domain. As a result, the surface shape of the measured object

*f*(

*x*,

*y*) is reconstructed by simply performing inverse Fourier transform to

*F*(

*u*,

*ν*).

## 3. Simulation results

To verify the viability of the proposed stereoscopic PMD for specular surface measurements, a simulation is carried out by measuring a square flat mirror.

For simulation, we need to settle the simulation model and to define the parameters of both the LCD and the cameras. For the LCD, we assume it is strictly homogeneous, and it has a resolution of 1000×1000 pixels with an LCD pixel size of 500×500 μm^{2}. What is more, to simplify the simulation process, we directly assume each LCD pixel has a same dimension. Nevertheless, we want to note that the real LCD size varies due to fabrication errors, and a pre-calibration is necessary if a high measurement accuracy is desired. On the other hand, two pinhole cameras are assumed for simulation disregarding any lens distortions. For both cameras, the CCD pixel size is set as 4.65×4.65 μm^{2}, the CCD resolution is 1000×1000 pixels, and the focal length of the employed lens is 12 mm. By setting these parameters, we let the simulation model be consistent to the real experiment set-up (see Section 4). Subsequently, we position the LCD directly facing the square flat mirror and we set the distance from the LCD to the mirror at 1950 mm. The stereoscopic PMD model for simulation is presented in Fig. 5. Both cameras are located 1800 mm from the flat mirror surface, and they are tilted to the surface at an angle of 3.8°. The two cameras are separated by 240 mm, with both camera pinholes and the LCD center located in a same plane (i.e., *XZ*-plane in Fig. 5).

Once the simulation model is determined, we assume the SFPs in the LCD can be reflected into any pinhole camera. Thus, we obtain the phase distributions in both pinhole camera image planes and the LCD plane. Next, we use phase difference minimization with second order polynomial fitting iteration to acquire the normal distribution with respect to the studied flat surface. Subsequently, we directly use 2D Fourier integration to reconstruct the surface, yielding the simulation result in Fig. 6(a). Simulated flat mirror in Fig. 6(a) has a dimension of 100×100 mm^{2}, and the surface demonstrates great flatness as the obtained peak-to-valley (PV) value is below 0.03 nm.

By considering the great surface homogeneity in Fig. 6(a), we prove that our proposed phase difference minimization algorithm is feasible to perform surface measurement in a stereoscopic deflectometry system free of calibration errors. Subsequently, we further study the influence of system calibration errors to the surface measurement results. First of all, let us consider the flat surface is slightly out of focus, where the distance between the LCD and the flat surface is 1905 mm rather than 1900 mm. Under this scenario, we use the same simulation model to measure the surface shape, and the result is given in Fig. 6(b). By comparing Fig. 6(b) to Fig. 6(a), we barely see any difference. Thus, we verify that even though the measured object is slightly out of focus, our stereoscopic PMD system still guarantees an accurate measurement. Afterwards, we evaluate the LCD pixel size error by setting its size to 500.01 μm instead of the ideal LCD pixel size of 500 μm. The simulation result with respect to LCD pixel error is given in Fig. 6(c), showing that the reconstructed surface deviates from a flat shape and the surface PV is deteriorated to nearly 10 μm. Therefore, a pre-calibration of LCD pixel size is required if a high accuracy is desired. In continuation, we also analyze the angular errors of LCD positioning. First, let us consider that the LCD is inclined with respect to *X* axis (see *O*-*X*,*Y*,*Z* in Fig. 5) by 0.1°, yielding the reconstructed surface shown in Fig. 6(d). Next, we also import an inclination error of 0.1° to both *Y* and *Z* axes, with the corresponding reconstructed surfaces shown in Figs. 6(e) and 6(f). It is obvious from the above simulations that inclination errors dramatically influence the surface reconstruction accuracy. It is shown that inclination in *Y* axis or *Z* axis is especially noticeable, as both cases result in PV errors greater than 100 μm.

## 4. Experimental measurements

To verify the feasibility of our proposed stereoscopic phase measuring deflectometry system, we experimentally calibrated the system and measured a flat and a spherical mirror. Measurements included the repeatability tests as well.

#### 4.1 Stereoscopic PMD system calibration

We have demonstrated a stereoscopic PMD system consists of an LCD screen to display sinusoidal fringe patterns, two cameras to capture the reflective patterns and the test surface. Three terms should be considered for system calibration, such as LCD gamma nonlinearity calibration, camera parameter calibration and system geometric calibration.

The LCD screen (HP Elite Display E231, resolution of 1920×1080 pixels, pixel size of 0.265×0.265 mm^{2}) implemented in our system presented gamma illumination nonlinearity, so it deteriorated the linear gray level-intensity response and therefore degraded the accuracy of phase retrieval. A lookup table was built [25] to correct the LCD illumination nonlinearity. The correction outcome was tested by performing the gamma test and it was found that gray level-intensity response demonstrated a linear tendency between the LCD and both cameras after gamma calibration.

Once the LCD gamma nonlinearity is calibrated, the camera calibration is subsequently performed. In our stereoscopic PMD system, we used the CCD cameras (charged-coupled device) provided by Basler (Basler Scout scA1000-30gc) with a resolution of 1034×779 pixels and a pixel size of 4.65×4.65 µm^{2}. The fixed focal length objective lenses were purchased from Thorlabs (MVL12M23) with a focal length of 12 mm. Each camera was provided with an individual adjustable support which contains holders for fiducial markers. Subsequently, both cameras were fixed on a breadboard facing the same direction and they were separated by nearly 240 mm. A well-known camera calibration toolbox [19,26] was used to perform the camera calibration, thus yielding camera intrinsic parameters such as focal length, principle point and distortion coefficients. Apart from that, we also obtained the geometric relation between the stereoscopic cameras by extrinsic parameter calibration.

Finally, we calibrated the geometric relation between the LCD and the cameras. In the first step, we located the breadboard, which was carrying both cameras, at position *A* (see Fig. 7). The optical center of Camera 1 was regarded as the coordinate system origin (see *O _{C}*-

*X*,

_{C}*Y*,Z

_{C}*in Fig. 7). Then, we attached three fiducial markers on each camera, and their position was determined using a laser tracker (FARON Vantage Laser Tracker) while cameras were located at position*

_{C}*A*. Subsequently, the breadboard with cameras was moved to position

*B*, and the positions of fiducial markers were measured again. In this way, the geometric transformation

**, of breadboard translation from**

*RT*_{BtA}*B*to

*A*is obtained, which is a 4×4 matrix.

We also generated a digital checkerboard pattern on the LCD panel in such way that the origin of the checkerboard coincided with the LCD coordinate system origin (see *O _{L}*-

*X*,

_{L}*Y*,

_{L}*Z*in Fig. 7). So, we used Camera 1 combined with camera calibration toolbox [26] to determine a 4×4 geometric transformation matrix from the LCD to Camera 1 at position

_{L}*B*as

**(see the green dash line in Fig. 7). Finally, we need to highlight that Camera 1 coordinate system was different from the breadboard coordinate system (this is not shown in Fig. 7). The transformation matrix from the breadboard coordinate system to Camera 1 coordinate system**

*RT*_{LtB}**(4×4 matrix) was established through hand-eye calibration [27]. Once we obtained all three geometric transformation relations, the coordinate system transformation from the LCD to Camera 1 were established through**

*HtE***, where**

*RT***=**

*RT***·**

*HtE***·**

*RT*_{BtA}**·**

*HtE*^{−1}**. Once this complete coordinate system transformation was determined, we were able to perform inverse ray tracing for surface normal determination.**

*RT*_{LtB}#### 4.2 Specular surface measurements

After the stereoscopic PMD system was calibrated, we implemented the specular surface measurements by measuring both a flat specular mirror and a spherical mirror.

For flat mirror measurement, we selected a broadband dielectric mirror of diameter 50.8 mm (Thorlabs, BB1-E02), and the mirror flatness of *λ*/10 at 633 nm. The mirror was located nearly 930 mm from the LCD screen. Mirror position was adjusted so that sinusoidal fringe patterns (SFPs) were reflected into both cameras. The experimental set-up is shown in Fig. 8. Subsequently, we generated on the LCD screen orthogonal SFPs with three different frequencies. The horizontal SFP frequencies were 540 pixels, 120 pixels and 30 pixels, and the vertical SFP frequencies were 480 pixels, 120 pixels and 30 pixels. Within each frequency, we introduced four patterns with a phase-shifting value of *π/*2. Once the SFPs reflected from the flat mirror were captured (see right side sub-figure in Fig. 8), phase unwrapping [22] was used to obtain the unwrapped phase maps in each camera. Finally, second order polynomial fitting phase difference minimization (as presented in Section 2.2) and 2D Fourier integration based surface reconstruction (as discussed in Section 2.3) were used to reconstruct the surface. Here, the phase difference minimization criteria *ε* was set to a value of 0.01. The resulting reconstructed surface shape of the flat mirror is shown in Fig. 9(a).

The reconstructed surface presented in Fig. 9(a) shows a smooth flatness, the surface root mean square (RMS) is 3.90 µm, and its peak-to-valley (PV) value is 11.54 µm. Upon further examination of the reconstructed surface, we find that the central section of mirror presents great flatness whereas a slight bending is observed in the peripheral surface areas close to edges. Note that such bending is in fact introduced by system calibration errors, which has been suggested by the simulation results in Fig. 6.

To determine the measurement accuracy, we also used an interferometer (Zygo Corp., MiniFiz 100) to measure this same flat mirror and the result is given in Fig. 9(b). The flat mirror measured by the interferometer presents a PV of 263.20 nm and an RMS of 18.58 nm. While the higher precision of interferometric measurement is obvious, we want to highlight again that the stereoscopic PMD measurement accuracy is mainly limited by the calibration errors. In fact, if we could use an LCD screen with better performance and cameras with negligible distortion, as well as we perform a more precise system calibration, the surface measurement accuracy could be further enhanced.

Besides the measurement accuracy, the measurement robustness is critical as well. Therefore, repeatability tests were performed with this same flat mirror to examine the system robustness. Specifically, the flat mirror was located at a same position and we measured it for four times using the same SFPs. Once we obtain these four measurement results, their average surface shape is calculated, and this average surface presented an RMS of 4.09 µm [see Fig. 10(a)]. Moreover, we also compared this average surface to any measurement results, and the surface shape difference maps are given in Fig. 10(b). The RMS of each sub-figure in Fig. 10(b) were as well calculated to represent the system repeatability test performance. Finally, we can easily notice that sub-figures in Fig. 10(b) only presented a maximum RMS of 35.10 nm, so it is obvious that stereoscopic PMD system has great measurement robustness.

The flat mirror measurement has justified the feasibility and robustness of stereoscopic PMD. We further measured a spherical mirror to verify that this system is also able to measure curved surfaces. Here, the spherical mirror was a dielectric-coated spherical concave mirror (Thorlabs, CM750-500-E02). The mirror has a curvature radius of 1000 mm with a diameter of 75 mm. The surface irregularity is *λ/*4 at 633 nm. Here, we want to emphasize that as this spherical mirror has a great curvature radius, it is difficult to be measured by an interferometer unless we use a unique and expensive large curvature reference mirror. We located the spherical mirror nearly 930 mm to the LCD screen, and we set the frequencies of both the horizontal and vertical SFP as 120 pixels, 30 pixels and 24 pixels. By using the same surface reconstruction method, we reconstruct the spherical mirror surface in Fig. 11(a), and it demonstrates great smoothness. What is more, we further select a horizontal cross line passing through the mirror center [from the middle of the mirror left edge to the middle of the mirror right edge, see red dash line in Fig. 11(a)], and we calculate its curvature radius as 1001.40 mm. Now, by comparing the experimentally measured curvature radius to an ideal spherical curve of radius 1000 mm, we find they are nearly identical as shown in Fig. 11(b). Therefore, the accuracy of phase difference minimization based stereoscopic PMD for curved surface measurement is verified. Finally, we compare the surface shape in Fig. 11(a) to an ideal spherical surface of a curvature radius of 1000 mm, and the surface difference map is given in Fig. 11(c). From Fig. 11(c), we notice that the measured surface is bended at the edge, where such characteristic also appears in flat mirror measurement [see Fig. 9(a)]. So, we can also attribute such surface bending to system calibration error.

Apart from the spherical surface shape measurement, repeatability tests were conducted as well. Similar to the flat mirror case, we measured the spherical mirror located at a same position using the same SFPs for four times. Once the four surface shapes were obtained, we calculated their average surface shape, which is shown in Fig. 12(a). Afterwards, we used this average surface map to calculate the difference to any distinct measurements. The resulting surface difference maps are presented in Fig. 12(b). The maximum RMS within all surface difference maps is 1.43 µm. Hence, once again the measurement robustness of stereoscopic PMD method is confirmed.

Finally, another repeatability test was carried out at different SFP frequencies. In this way, we can assess how different SFP frequencies will influence the measurement robustness. Frequencies of both the horizontal and vertical SFPs were changed to 120 pixels, 30 pixels and 20 pixels (the previous measurement used the frequencies of 120 pixels, 30 pixels and 24 pixels), and the repeatability test was carried out again. Results can be seen in Fig. 13. Note that both the averaged surface shape and the maximum surface difference map RMS are nearly identical to the results shown in Fig. 12, and the maximum RMS is only 2.49 *µ**m*. Thus, by using different SFP frequencies we further corroborate the measurement robustness of our proposed stereoscopic PMD method.

## 5. Conclusions

In this paper, we demonstrate a phase difference minimization algorithm based stereoscopic phase measuring deflectometry (PMD) system to measure the surface shape of a specular object. By introducing stereo cameras, the height-normal ambiguity is solved without introducing any system component movements once the system is calibrated. More importantly, the phase difference minimization method, which is executed by using the phase distributions between the LCD screen and the camera image plane, is proposed to determine the surface point. Moreover, for the first time, an efficient second order polynomial fitting method is suggested in this work for phase difference minimization. Afterwards, the surface shape is reconstructed by performing both bi-cubic sub-pixel interpolation and 2D Fourier integration. To verify the feasibility of our stereoscopic PMD system, a flat mirror and a spherical mirror were measured. The flat mirror measurement yields surface RMS value of only 3.90 μm and measurement repeatability of 35.10 nm, thus showing high measurement precision and great measurement robustness. The spherical mirror measurement, on the other hand, presents a measured curvature radius of 1001.4 mm, which is identical to 1000 mm specified by Thorlabs. Apart from that, measurement repeatability was also tested using two different series of SFP frequencies. Both tests resulted in maximum repeatability RMS values of 1.43 µm and 2.49 µm, respectively. Experimental measurement results show that we have comprehensively verified the feasibility and efficiency of our proposed phase difference minimization based stereoscopic phase measuring deflectometry system. Future research will focus on system calibration error elimination and measurement speed enhancement.

## Funding

Agència de Gestió d'Ajuts Universitaris i de Recerca (SGR 2017-001500); Ministerio de Economía, Industria y Competitividad, Gobierno de España (Fondos Feder, RTI2018-097107-B-C31); Federación Española de Enfermedades Raras (Spanish MinCIU, RTI2018- 097107-B-C32).

## Acknowledgments

The authors thank Dr. Lei Huang from Brookhaven National Laboratory for his suggestions on camera calibration, Dr. Pablo Pedrera from ALBA Synchrotron Light Source for organizing a valuable laboratory environment. H. Zhang thanks Dr. Daniel Lau from University of Kentucky for his useful advice on gamma calibration.

## Disclosures

The authors declare no conflict of interest.

## References

**1. **D. Malacara, * Optical Shop Testing* (Wiley, 2007).

**2. **M. C. Knauer, J. Kaminski, and G. Hausler, “Phase measuring deflectometry: a new approach to measure specular free-form surfaces,” Proc. SPIE **5457**, 366 (2004). [CrossRef]

**3. **Y. Liu, S. Huang, Z. Zhang, N. Gao, F. Gao, and X. Jiang, “Full-field 3D shape measurement of discontinuous specular objects by direct phase measuring deflectometry,” Sci. Rep. **7**(1), 10293 (2017). [CrossRef]

**4. **Y. Tang, X. Su, Y. Liu, and H. Jing, “3D shape measurement of the aspheric mirror by advanced phase measuring deflectometry,” Opt. Express **16**(19), 15090–15096 (2008). [CrossRef]

**5. **Y. Tang, X. Su, F. Wu, and Y. Liu, “A novel phase measuring deflectometry for aspheric mirror test,” Opt. Express **17**(22), 19778–19784 (2009). [CrossRef]

**6. **P. Su, R. E. Parks, L. Wang, R. P. Angel, and J. H. Burge, “Software configurable optical test system: a computerized reverse Hartmann test,” Appl. Opt. **49**(23), 4404–4412 (2010). [CrossRef]

**7. **P. Su, Y. Wang, J. H. Burge, K. Kaznatcheev, and M. Idir, “Non-null full field X-ray mirror metrology using SCOTS: a reflection deflectometry approach,” Opt. Express **20**(11), 12393–12406 (2012). [CrossRef]

**8. **Z. Niu, X. Xu, X. Zhang, W. Wang, Y. Zhu, J. Ye, M. Xu, and X. Jiang, “Efficient phase retrieval of two-directional phase-shifting fringe patterns using geometric constraints of deflectometry,” Opt. Express **27**(6), 8195–8207 (2019). [CrossRef]

**9. **L. Huang, C. S. Ng, and A. K. Asundi, “Dynamic three-dimensional sensing for specular surface with monoscopic fringe reflectometry,” Opt. Express **19**(13), 12809–12814 (2011). [CrossRef]

**10. **T. Bothe, W. Li, C. von Kopylow, and W. P. O. Juptner, “High-resolution 3D shape measurement on specular surfaces by fringe reflection,” Proc. SPIE **5457**, 411 (2004). [CrossRef]

**11. **X. Xu, X. Zhang, Z. Niu, W. Wang, Y. Zhu, and M. Xu, “Self-calibration of in situ monoscopic deflectometric measurement in precision optical manufacturing,” Opt. Express **27**(5), 7523–7536 (2019). [CrossRef]

**12. **L. Huang, J. Xue, B. Gao, C. McPherson, J. Beverage, and M. Idir, “Modal phase measuring deflectometry,” Opt. Express **24**(21), 24649–24664 (2016). [CrossRef]

**13. **Z. Niu, N. Gao, Z. Zhang, F. Gao, and X. Jiang, “3D shape measurement of discontinuous specular objects based on advanced PMD with bi-telecentric lens,” Opt. Express **26**(2), 1615–1632 (2018). [CrossRef]

**14. **C. Li, Y. Li, X. Xiao, X. Zhang, and D. Tu, “Phase measurement deflectometry with refraction model and its calibration,” Opt. Express **26**(26), 33510–33522 (2018). [CrossRef]

**15. **T. Tao, Q. Chen, J. Da, S. Feng, Y. Hu, and C. Zuo, “Real-time 3-D shape measurement with composite phase-shifting fringes and multi-view system,” Opt. Express **24**(18), 20253–20269 (2016). [CrossRef]

**16. **S. Heist, P. Dietrich, M. Landmann, P. Kühmstedt, G. Notni, and A. Tünnermann, “GOBO projection for 3D measurement at highest frame rates: a performance analysis,” Light: Sci. Appl. **7**(1), 71 (2018). [CrossRef]

**17. **L. P. Yaroslavsky, A. Moreno, and J. Campos, “Frequency responses and resolving power of numerical integration of sampled data,” Opt. Express **13**(8), 2892–2905 (2005). [CrossRef]

**18. **C. Zuo, S. Feng, L. Huang, T. Tao, W. Yin, and Q. Chen, “Phase shifting algorithms for fringe projection profilometry: A review,” Opt. Laser Eng. **109**, 23–59 (2018). [CrossRef]

**19. **J. Heikkila and O. Silven, “A four-step camera calibration procedure with implicit image correction,” in Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition (1997), pp. 1106–1112.

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

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

**22. **C. Zuo, L. Huang, M. Zhang, Q. Chen, and A. Asundi, “Temporal phase unwrapping algorithms for fringe projection profilometry: A comparative review,” Opt. Laser Eng. **85**, 84–103 (2016). [CrossRef]

**23. **B. Zhao and Y. Surrel, “Effect of quantization error on the computed phase of phase-shifting measurements,” Appl. Opt. **36**(10), 2070–2075 (1997). [CrossRef]

**24. **W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, * Numerical recipes in C: the art of scientific computing* (Cambridge University, 1992), Chap. 3.

**25. **P. S. Huang, C. Zhang, and F.-P. Chiang, “High-speed 3-D shape measurement based on digital fringe projection,” Opt. Eng. **42**(1), 163–169 (2003). [CrossRef]

**26. **J. Y. Bouguet, “Camera calibration toolbox for matlab,” http://www.vision.caltech.edu/bouguetj/calib_doc.

**27. **R. Y. Tsai and R. K. Lenz, “A new technique for fully autonomous and efficient 3D robotics hand/eye calibration,” IEEE Trans. Robot. Automat. **5**(3), 345–358 (1989). [CrossRef]