An Improved Algorithm for 3D Reconstruction Integration Based on Stripe Reflection Method

Xu Gao, Peiji Guo^*

School of Optoelectronic Science and Engineering, Key Lab of Advanced Optical Manufacturing Technologies of Jiangsu Province and Key Lab of Modern Optical Technologies of Education Ministry of China, Soochow University, Suzhou, China.

DOI: 10.4236/opj.2019.98B002                

Abstract

This paper introduces the basic principle of stripe reflection method and proposes an improved algorithm on the traditional Southwell gradient iterative integration algorithm. The algorithm adds a coefficient value with an attenuation factor to the compensation height value and the value of the attenuation factor is changed by the determination of the compensation height threshold. Through computer simulation, the fitting error of the reconstructed surface show that the RMS of the new method is one order of magnitude better than the traditional algorithm and the PV value of the high frequency part is about 1/15 of the traditional algorithm. It is proved that the improved algorithm can effectively improve the convergence and noise resistance of the iterative algorithm.

Keywords

Stripe Reflection Method, Gradient Iterative Integration, Fitting Error

1. Introduction

The fringe reflection method [1] is based on the simple principle of light reflection. The method has the advantages of simple structure, high sensitivity, large dynamic range, and is generally used for full-caliber testing of the aspherical surface shape during the rough polishing stage. Many experts and scholars have carried out many researches and practices on this technology. Su Peng et al. of the University of Arizona in the United States proposed a SCOTS testing method and successfully detected the 8.2 m GMT sub-mirror through high-precision calibration [2] [3]. The Sichuan University [4] [5] and the Chinese Academy of Sciences [6] have done a lot of research on the stripe reflection technology and have achieved certain results. This paper introduces the basic principle of the stripe reflection method and introduces an iterative algorithm based on the Southwell integration method. Then we propose a compensation algorithm with attenuation factor parameters based on the traditional iterative integration. The computer simulation results show that the new method has better convergence in a single iteration compared with the traditional method. The algorithm for changing the coefficient has a better fitting result than the algorithm without changing the coefficient, and verifies the superiority of the improved algorithm.

2. The Basic Principle of the Stripe Reflection Method

The measurement model is shown in Figure 1. The computer produces a standard sinusoidal stripe that transmits the encoded fringe image onto the LCD screen. The stripe displayed on the screen is projected onto the reflective surface. The deformed fringe image is captured by a CCD camera, which records the gradient and height of each reflective point on the reflective surface.

The intensity of each point on the deformed fringe pattern after reflection by the mirror to be tested can usually be expressed as:

$I_n(x,y)=a(x,y)+b(x,y)\cos [\phi (x,y)+{\sigma }_n]$ (1)

where $I_n\left(x,y\right)$ is the light intensity of each point of the stripe diagram, $a\left(x,y\right)$ is the background light intensity, $b(x,y)$ is amplitude modulation, and $\phi \left(x,y\right)$ is the phase modulation to be solved, and ${\sigma }_n$ is the additional phase, which is determined when the fringes are generated. The computer generated grating can accurately control the magnitude of the single phase shift. The calculation of the phase modulation $\phi \left(x,y\right)$ is shown in Equation (2)

$\phi (x,y)=-{\tan }^{-1}\left[\frac{{\displaystyle \underset{n-1}{\overset{N}{\sum }}{I}_n(x,y)\sin {\sigma }_n}}{{\displaystyle \underset{n-1}{\overset{N}{\sum }}{I}_n(x,y)\cos {\sigma }_n}}\right],n=1,2,\mathrm{...},N$ (2)

The arctangent function obtains the truncated phase of [−π, π). In order to determine the one-to-one correspondence between the pixel of the screen and the pixel of the mirror, it is necessary to expand the phase to form a continuous phase. The method of phase unwrapping includes the spatial phase expansion method [7] [8] and time phase expansion method [9] [10]. The former mainly

Figure 1. Stripe reflection measurement model.

includes branch cutting method, quality map guiding method, minimum norm method; the latter mainly includes binary encoding and multi-frequency heterodyne method.

3. Gradient Integral Algorithm Analysis

Measuring the surface of the object by the stripe reflection system, we can obtain two gradient distributions in the mutually perpendicular direction. The gradient refers to the partial derivatives $g_x(x,y)$ and $g_y(x,y)$ obtained by the height $z(x,y)$ of the surface in the directions of x and y. The calculation is as shown in Equation (3) and Equation (4):

$g_x(x,y)=\frac{\partial z(x,y)}{\partial x}$ , (3)

$g_y(x,y)=\frac{\partial z(x,y)}{\partial y}$ , (4)

In order to obtain the specular surface 3D height information, we need to integrate the derivative:

$z(x,y)={\displaystyle \int g_{x}dx+{\displaystyle \int g_{y}dy}}$ , (5)

In the actual testing process, the gradient field is not a conservative field due to the influence of system error and image noise. Southwell model based on regional wave front reconstruction method can not only effectively suppress noise but also deal with complex connected regions. Therefore, it has become the mainstream algorithm of gradient integral.

3.1. Southwell Model Algorithm

The Southwell grid model is shown in Figure 2. The matrix size is M × N. The solid black dot in the graph represents the height point and the arrow represents the gradient direction of the point x, y. The average of two adjacent gradient values is equal to the ratio of the difference in height to the pixel width (h).

The gradient matrix is the same size as the height matrix. The relationship between the gradient data and the height data is only related to the same row or the same column. The above relationship is expressed as

$\{\begin{array}{l}\frac{1}{2}\left(g_{i,j+1}^x+g_{i,j}^x\right)=\frac{1}{h}(z_{i,j+1}-z_{i,j})i=1,2,\mathrm{...},M;j=1,2,\mathrm{...},N-1\\ \frac{1}{2}\left(g_{i+1,j}^y+g_{i,j}^y\right)=\frac{1}{h}(z_{i+1,j}-z_{i,j})i=1,2,\mathrm{...},M-1;j=1,2,\mathrm{...},N\end{array}$ (6)

Z represents a height matrix, the size is MN × 1; D represents a coefficient matrix, which is a sparse matrix size [(M − 1)N + M(N − 1)] × MN; G represents a gradient matrix, and the size [(M − 1)N + M(N − 1)] × 1. Equation (6) can be rewritten into:

$DZ=G$ , (7)

Generally, the surface matrix is relatively large. In this case, the D matrix is a singular matrix. The augmented matrix $D_f$ needs to be replaced by the matrix

Figure 2. Southwell grid model.

D, and the least norm solution of Z can be obtained by using the least square method:

$Z={\left(D_f^T{D}_f\right)}^{-1}{D}_f^{T}G$ . (8)

3.2. Improved Iterative Algorithm Based on Southwell Model

The rate of convergence and the ability to suppress noise in the late iteration are two conditions for evaluating the merits of an iterative algorithm. Especially for large-caliber mirrors, the gradient data needs to be reprocessed for each iteration. The higher the convergence efficiency, the fewer the number of iterations. As the number of iterations increases, the ability to suppress gradient noise in the iteration ensures the accuracy of the final fitted shape.

The traditional iterative algorithm is showed as follows [11]

1) Using the traditional integration method to obtain the original surface height $z^0(x,y)$ as the initial value of the iteration;

2) Calculate the derivatives $q(x,y)$ and $p(x,y)$ of the current height $z^n(x,y)$ in the x and y directions, and make a difference with the original gradient to obtain gradient residuals $d{g}^x(x,y)$ and $d{g}^y(x,y)$ ;

3) Using traditional integration method to integrate $d{g}^x(x,y)$ and $d{g}^y(x,y)$ to obtain the compensation height value $z_c^n(x,y)$ ;

4) The compensated height values $z^n\left(x,y\right)=z^{n-1}\left(x,y\right)+z_c^n(x,y)/n_k$ , n_k = 3, 4.0909, 4.9476, 5.6768 (k = 1, 2, 3, 4);

5) Repeat the iterative process of 2-4 until the iterations of adjacent two iterations are less than the threshold, $\max \left|z^n-z^{n-1}\right|$ < threshold.

This section presents an improved algorithm that includes two improvements. First, use the two variables of attenuation factor t and iteration number n to control the contribution value of the compensation height to participate in the iteration. Second, the threshold is set for the nth compensation height during iteration. The formula of the new iterative algorithm can be written as:

$z^n(x,y)=z^{n-1}(x,y)+\omega \cdot z_c^n(x,y)$ , (9)

$\omega =\frac{1}{{(n-1)}^t}$ , (10)

The height of the nth compensation is related to the height threshold:

$\text{RMS}\left(z_c^n\right)<\frac{\text{threshold}}{n\cdot \text{m}}$ , (11)

The value of m is generally between [1] [5], Coefficient $\omega$ determines the weight of the compensation height value to participate in the iteration. In the initial stage of the iteration, a smaller $t$ can increase the convergence rate. With the progress of iteration, in order to reduce excessive noise participation in the iteration, we set the RMS threshold of the nth compensation height. If it is larger than the threshold value, the value of attenuation factor $t$ is increased to prevent excessive noise data in the compensation height from affecting the fitting results. Through computer simulation, the traditional algorithm and the improved algorithm iteratively integrate the gradient data obtained by the stripe reflection measurement system, and analyze the wavefront reconstruction accuracy.

4. Computer Simulation

The simulated face is divided into two parts: $z_1$ is the standard paraboloid:

$z_1=\frac{x^2}{50}+\frac{y^2}{50}$ , (12)

$z_2$ is set to a high-order shape in a certain area:

$z_2=\frac{1}{0.5{\cos }^2(x^6+y^4)+{(0.25{(0.5x-2)}^4+2)}^2+0.125{({(0.5y)}^5+2)}^2}$ , (13)

The size of $z_1$ and $z_2$ is 501 × 501 pixels, and the pixel size is 0.04 mm. In order to get closer to the actual testing conditions, we combined the gradient data of surface shape with the gaussian noise whose SNR = 100:1. Using the traditional Southwell integral $z_0$ as the initial shape, after four iterations we compare the high frequency partial fit results on the green lines in Figure 3. The simulation results are shown in Figure 4.

In order to compare the capability of the three iterative methods to iterative errors in the high-frequency part, we gather the high-frequency errors fitted by the three methods into the same figure, as shown in Figure 5, and list the RMS of the four iterations, as shown in Table 1.

In order to verify the effect of controlling the attenuation factor by setting compensation high threshold in the later iteration, we will perform 10 iterations without changing the initial attenuation coefficient t and controlling the attenuation coefficient, and compare the rms results of the last four fittings. The result is shown in Figure 6. The height average of the last four compensations is shown in Table 2.

5. Summary and Discussion

This paper introduces the basic principle of the stripe reflection method, and proposes an improved algorithm based on traditional Southwell iterative

(a) (b)

Figure 3. The simulation surface shape. (a) The two-dimensional shape of $z_2$ , (b) Two-dimensional shape synthesized by $z_1$ and $z_2$ .

(a) (b) (c) (d) (e) (f)

Figure 4. Gradient integral method reconstruction error and high frequency fitting error. (a) Traditional Southwell integral surface residual error, (b) Traditional Southwell integral high frequency line fitting error, (c) Traditional iterative integral surface residual error, (d) Traditional iterative integration high frequency line fitting error, (e) Improved surface residual error of iterative integration, (f) Improved iterative integral high frequency line fitting error.

Figure 5. Comparison of high-frequency line error distribution.

Figure 6. Root mean square error comparison.

integration. The algorithm uses the attenuation coefficient to control the compensation value and the attenuation coefficient t is controlled by the compensation height threshold. Through computer simulation, it is verified that the algorithm has better convergence rate in the early iteration period and better fitting

Table 1. The error size of single iteration reconstruction in iteration process.

Table 2. The average compensation height of the last four iterations (mm).

of high frequency detail. In the latter part of the iteration, by changing the size of t, it can effectively prevent too many noise points from participating in the iteration and leading to poor surface performance. The simulation results show that the improved algorithm effectively suppresses the influence of noise and has better fitting accuracy than the traditional algorithm at low SNR.

Funding

Supported by the NSFC (NO.61378055), Advantages of disciplines in colleges and universities in Jiangsu Province construction grant program.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Hung, Y.Y., Lin, L., Shang, H.M., et al. (2000) Practical Three-Dimensional Computer Vision Techniques for Full-Field Surface Measurement. Opt. Eng., 39, 143-149. https://doi.org/10.1117/1.602345
[2]	Su, P., Parks, R.E., Wang, L.R., et al. (2010) Software Configurable Optical Test System: A Computerized Reverse Hartmann Test. Applied Optics, 49, 4404–4412. https://doi.org/10.1364/ao.49.004404
[3]	Huang, R., Su, P. and Burgea, J.H. (2015) Deflectometry Measurement of Daniel K. Inouye Solar Telescope Primary Mirror. Proceedings of the SPIE, 9575, Article ID: 957515.
[4]	Zhao, W.C., Fan, B., Wu, F., et al. (2013) Experimental Analysis of Reflector Test Based on Phase Measuring Deflectometry. Acta Optica Sinica, 33, Article ID: 0112002. https://doi.org/10.3788/aos201333.0112002
[5]	Liu, Y.K., Su, X.Y. and Wu, Q.Y. (2006) Three-Dimensional Shape Measurement for Specular Surface Based on Fringe Reflection. Acta Optica Sinica, 26, 1636-1640.
[6]	Wang, H.R., Li, B., Wang, Z.F., et al. (2013) Plane Measurement of Trough Parabolic Element Mirror Based on Stripe Reflection Method. Acta Optica Sinica, 33, Article ID: 0112007.
[7]	Strand, J., Taxt, T. and Jain, A.K. (1999) Two-Dimensional Phase Unwrapping Using a Block Least-Square Method. IEEE Tran. Image Process, 8, 375-386. https://doi.org/10.1109/83.748892
[8]	Goldstein, R.M., Zevker, H.A. and Werner, C.L. (1988) Statellite Radar Interferometry: Two-Dimensional Phase Unwrapping. Radio Science, 23, 713-720. https://doi.org/10.1029/rs023i004p00713
[9]	Strand, J., Taxt, T. and Jain. A.K. (1999) Two-Dimensional Phase Unwrapping Using a Block Least-Square Method. IEEE Tran. Image Process, 8: 375-386. https://doi.org/10.1109/83.748892
[10]	Huntley, J.M. and Saldner, H.O. (1993) Temporal Phase-Unwrapping Algorithm for Automated Interferogram Analysis. Applied Optics, 32, 3047-3052. https://doi.org/10.1364/ao.32.003047
[11]	Huang, L. and Asundi, A. (2012) Improvement of Least-Squares Integration Method with Iterative Compensations in Fringe Reflectometry. Applied Optics, 51, 7459-7465. https://doi.org/10.1364/ao.51.007459