Improving Least-Squares Surface Reconstruction through Fourth-Order Adams Method and Iterative Compensation

Geng An¹, Fenlin Yang^2*, Guangying Liu², Feifan Fu²

¹School of Communication and Electronic Engineering, Jishou University, Jishou, China.
²School of Mathematics and Statistics, Jishou University, Jishou, China.

DOI: 10.4236/oalib.1110454                

Abstract

The Southwell model stands as a prominent algorithm within the realm of the least squares surface reconstruction, finding wide application. This algorithm boasts notable merits, including rapid computation and an approximation of the reconstructed surface that closely mirrors reality. Nevertheless, it is not without its drawbacks, as it exhibits substantial reconstruction errors and proves to be susceptible to the presence of noisy data. To enhance the precision of the reconstructed object’s three-dimensional surface, this paper puts forth an enhanced least squares surface reconstruction algorithm based on the fourth-order Adams method and iterative compensation. Initially, the fourth-order Adams method is employed to establish the connection between the measured gradient and the unknown surface height in the Southwell model. Subsequently, Tikhonov regularization is introduced to mitigate the impact of noise on the model. Ultimately, the accuracy is augmented through the utilization of an iterative compensation technique. Simulation experiments substantiate that, in comparison to alternative Southwell model algorithms, the proposed algorithm exhibits reduced time consumption and superior surface fitting accuracy.

Keywords

Least Square, Southwell Model, Adams, Surface Reconstruction

1. Introduction

The utilization of gradient integration algorithms enables the reconstruction of object surface morphology through the application of discrete gradient data. This approach has garnered considerable prominence in various domains, including medical imaging, 3D reconstruction, and facial recognition in recent years [1] [2] [3] . Currently, gradient integration algorithms can be categorized into three main approaches: the cross-path method, the Fourier transform integration method, and the regional wavefront reconstruction method. The cross-path method involves integrating the measured gradient data increments along both the X and Y axes. In contrast, the Fourier transform integration method adopts a global-based integration strategy, establishing the relationship between height data and gradient data in the frequency domain and obtaining height data through the Fourier inverse transform. On the other hand, the regional wavefront reconstruction method associates the measured gradient values with the unknown wavefront height values, ultimately estimating the surface height using the least squares method.

Based on the correspondence between gradient data and height values for each point on the reconstructed surface, the regional wavefront reconstruction method encompasses the Hudgin model [4] , the Fried model [5] , and the Southwell model [6] . Specifically, the Southwell model exhibits a gradient matrix that aligns in size with the height matrix, mirroring real measurement conditions. Notably, this model surpasses the other two counterparts in terms of measurement and reconstruction accuracy while maintaining computational simplicity. Consequently, its widespread adoption is justified.

To enhance the precision of the Southwell model, Li et al. [7] employed the one-dimensional Taylor theory to diminish truncation errors by augmenting the gradient point count within the operator. Ren et al. [8] proposed a method to reconstruct the surface by combining least squares integration and radial basis function integration under the assumption that the normal vector is perpendicular to the vector determined by the two endpoints. Huang et al. [9] introduced a more generalized difference operator by considering the relationship between the combination of gradients in the X and Y directions and the adjacent height difference.

Based on the preceding literature analysis, it becomes evident that various approaches have been employed to reconfigure the association between gradient data and height data within the Southwell model, alongside augmenting its constraints. However, these endeavors have fallen short in adequately addressing the impact of noise on the model, as well as rectifying the presence of errors. Consequently, the established model exhibits inherent limitations. To address these issues, we propose utilizing the fourth-order Adams method to establish a more robust connection between the measured gradient data and the unknown height values. Additionally, we introduce Tikhonov regularization as a means to mitigate the impact of noise. Furthermore, to enhance the overall reconstruction accuracy, we incorporate an iterative compensation technique.

2. Southwell Model

The Southwell model postulates an equivalence between the measured gradient position and the corresponding height position. Figure 1 illustrates the model’s schematic diagram, depicting the vertical distribution of measured gradients represented by arrows, while the dots signify the positions of the height values to be reconstructed.

The Southwell model achieves the third-order precision reconstruction through finite difference equations, which can be expressed as:

$\{\begin{array}{ll}{z}_{i,j+1}-z_{i,j}=\frac{1}{2}\left(g_{i,j+1}^x+g_{i,j}^x\right)\left(x_{i,j+1}-x_{i,j}\right)\hfill & i=1,2,\cdots ,M;j=1,2,\cdots ,N-1.\hfill \\ z_{i+1,j}-z_{i,j}=\frac{1}{2}\left(g_{i+1,j}^y+g_{i,j}^y\right)\left(y_{i+1,j}-y_{i,j}\right)\hfill & i=1,2,\cdots ,M-1;j=1,2,\cdots ,N.\hfill \end{array}$ (1)

In Equation (1): z indicates the height of the reconstruction point; $g^x\mathrm{,}{g}^y$ is the gradient value calculated in the corresponding direction; $x\mathrm{,}y$ is the coordinate in the corresponding direction; $M\mathrm{,}N$ indicates the length and width of the surface to be built. For the convenience of calculation, Equation (1) is written in matrix form as:

$DZ=G,$ (2)

where D is the coefficient matrix, which is a sparse matrix of size $\left[\left(M-1\right)N+M\left(N-1\right)\right]\times MN$ , which can be expressed as:

$D=\left[\begin{array}{ccccccccc}-1& 0& \cdots & 0& 1& 0& \cdots & \cdots & 0\\ 0& -1& 0& \cdots & 0& 1& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& \cdots & \cdots & 0& -1& 0& \cdots & 0& 1\\ -1& 1& 0& \cdots & \cdots & \cdots & \cdots & \cdots & 0\\ 0& -1& 1& 0& \cdots & \cdots & \cdots & \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& \cdots & \cdots & \cdots & \cdots & \cdots & 0& -1& 1\end{array}\right],$ (3)

Figure 1. Wavefront reconstructionr model for Southwell region.

Z is an unknown height matrix of size $MN\times 1$ , which can be expressed as:

$Z=\left[\begin{array}{c}{z}_{\mathrm{1,1}}\\ z_{\mathrm{2,1}}\\ ⋮\\ z_{M-\mathrm{1,}N}\\ z_{M\mathrm{,}N}\end{array}\right]\mathrm{,}$ (4)

G is a gradient matrix composed of measured gradient data in the x and y directions, with a size of $\left[\left(M-1\right)N+M\left(N-1\right)\right]\times 1$ , which can be expressed as:

$G=\frac{1}{2}\left[\begin{array}{c}\left(g_{\mathrm{1,2}}^x+g_{\mathrm{1,1}}^x\right)\left(x_{\mathrm{1,2}}-x_{\mathrm{1,1}}\right)\\ \left(g_{\mathrm{2,2}}^x+g_{\mathrm{2,1}}^x\right)\left(x_{\mathrm{2,2}}-x_{\mathrm{2,1}}\right)\\ ⋮\\ \left(g_{M\mathrm{,}N}^x+g_{M\mathrm{,}N-1}^x\right)\left(x_{M\mathrm{,}N}-x_{M\mathrm{,}N-1}\right)\\ \left(g_{\mathrm{2,1}}^y+g_{\mathrm{1,1}}^y\right)\left(y_{\mathrm{2,1}}-y_{\mathrm{1,1}}\right)\\ \left(g_{\mathrm{3,1}}^y+g_{\mathrm{2,1}}^y\right)\left(y_{\mathrm{3,1}}-y_{\mathrm{2,1}}\right)\\ ⋮\\ \left(g_{M\mathrm{,}N}^y+g_{M-\mathrm{1,}N}^y\right)\left(y_{M\mathrm{,}N}-y_{M-\mathrm{1,}N}\right)\end{array}\right]\mathrm{,}$ (5)

Since the number of equations is $\left(M-1\right)N+M\left(N-1\right)$ , and the number of height points to be solved is MN, the number of equations is greater than the number of unknowns, which belongs to over-determined equations, and the final least squares solution of the equation is:

$\left(D^{\text{T}}D\right)Z=D^{\text{T}}G.$ (6)

The gradient data within the Southwell model displays a localized correlation with the height data, subsequently undergoing global optimization through least squares. Nevertheless, the Southwell model is subject to inherent limitations, notably in its assumption of a quadratic linear relationship between height and gradient. Consequently, a single integration yields solely a quadratic surface shape, lacking the capacity to accurately fit detailed surface information. Furthermore, the Southwell model disperses the impact of noise across the entire reconstruction, yet fails to entirely eradicate its influence. In the subsequent section, we shall address and enhance the deficiencies inherent in the Southwell model.

3. Improved Model

3.1. Fourth Order Adams

The reconstruction accuracy of the renowned Southwell model is constrained by the quadratic term governing the height variation between two sampling points. In this study, we employ the fourth-order Adams method to reconstruct the association between the measured gradient data and the unknown height. Compared to the conventional Southwell model, the fourth-order Adams method exhibits a reduced truncation error, yielding superior performance.

Furthermore, the reconstruction outcomes are influenced by the constraints imposed on adjacent point heights. The classic Southwell model only considers the relationship between the current height and the height of the preceding point, overlooking constraints on neighboring points’ heights. Consequently, the relative positions of adjacent points become imprecise after reconstruction. To enhance the restoration of high-frequency surface details, we extend the adjacent point constraint to encompass three points. The fourth-order Adams method is outlined as follows:

$y_{n+1}=y_n+\frac{h}{24}\left(9f_{n+1}+19f_n-5f_{n-1}+f_{n-2}\right)\mathrm{,}$ (7)

Combined with the fourth-order Adams, the unknown height data can be expressed as:

$\{\begin{array}{c}{z}_{i,j+1}-z_{i,j}=\frac{9g_{i,j+1}^x+19g_{i,j}^x-5g_{i,j-1}^x+g_{i,j-2}^x}{24}\left(x_{i,j+1}-x_{i,j}\right)\\ i=1,2,\cdots ,M;j=3,4,\cdots ,N-1\\ z_{i+1,j}-z_{i,j}=\frac{9g_{i+1,j}^y+19g_{i,j}^y-5g_{i-1,j}^y+g_{i-2,j}^y}{24}\left(y_{i+1,j}-y_{i,j}\right)\\ i=3,4,\cdots ,M-1;j=1,2,\cdots ,N\end{array}$ (8)

At the same time, the above formula will lead to the problem of missing values at the boundary of the gradient matrix, and we use the traditional Southwell model to fill in the missing values. Finally, it is converted into a matrix form, which is the same as Equation (2), but the gradient matrix G should be expressed as:

$G=\frac{1}{24}\left[\begin{array}{c}12\left(g_{\mathrm{1,2}}^x+g_{\mathrm{1,1}}^x\right)\left(x_{\mathrm{1,2}}-x_{\mathrm{1,1}}\right)\\ 12\left(g_{\mathrm{2,2}}^x+g_{\mathrm{2,1}}^x\right)\left(x_{\mathrm{2,2}}-x_{\mathrm{2,1}}\right)\\ ⋮\\ \left(9g_{M\mathrm{,}N}^x+19g_{M\mathrm{,}N-1}^x-5g_{M\mathrm{,}N-2}^x+g_{M\mathrm{,}N-3}^x\right)\left(x_{M\mathrm{,}N}-x_{M\mathrm{,}N-1}\right)\\ 12\left(g_{\mathrm{2,1}}^y+g_{\mathrm{1,1}}^y\right)\left(y_{\mathrm{2,1}}-y_{\mathrm{1,1}}\right)\\ 12\left(g_{\mathrm{3,1}}^y+g_{\mathrm{2,1}}^y\right)\left(y_{\mathrm{3,1}}-y_{\mathrm{2,1}}\right)\\ ⋮\\ \left(9g_{M\mathrm{,}N}^y+19g_{M-\mathrm{1,}N}^y-5g_{M-\mathrm{2,}N}^y+g_{M-\mathrm{3,}N}^y\right)\left(y_{M\mathrm{,}N}-y_{M-\mathrm{1,}N}\right)\end{array}\right]\mathrm{.}$ (9)

3.2. Tikhonov Regularization

During the process of actual measurement, discrepancies between surface height and true height arise due to the influence of instrument error or noise affecting the measurement gradient. To express this phenomenon, we introduce the variable G, which represents the disturbed gradient data affected by errors and can be formulated as follows:

$G=G_{\text{tree}}+\eta \mathrm{,}$ (10)

Here, $G_{\text{true}}$ denotes the error-free gradient data, while $\eta$ represents the measured gradient error. In order to address the ill-conditioning issue associated with matrix D and the presence of measurement error $\eta$ in G, we introduce the generalized inverse matrix of D, denoted as $D^{†}$ . Consequently, the model solution can be expressed as $D^{†}G=D^{†}{G}_{\text{true}}+D^{†}\eta =Z_{\text{true}}+D^{†}\eta$ , taking into account the propagation error $D^{†}\eta$ . However, this model solution fails to accurately approximate $Z_{\text{true}}$ . To address this problem, one commonly employed approach is Tikhonov regularization, which allows us to obtain a minimization problem given by:

$\min \left\{{\Vert DZ-G\Vert }_2^2+\mu {\Vert Z\Vert }_2^2\right\}\mathrm{.}$ (11)

Here, $\mu >0$ serves as a regularization parameter, and the corresponding matrix equation is given by:

$\left(D^{\text{T}}D+\mu I\right)Z=D^{\text{T}}G.$ (12)

It can be observed that Equation (11) possesses a unique solution:

$Z_{\mu }={\left(D^{\text{T}}D+\mu I\right)}^{-1}{D}^{\text{T}}G.$ (13)

The Tikhonov minimization problem (11) consists of two terms: the first term, known as the fidelity term, ensures that the height approximates the observed surface height, while the second term, referred to as the penalty term, controls the smoothness of the reconstructed surface. The balance between these two terms is determined by the regularization parameter $\mu$ . If the choice of $\mu$ is inaccurate, $Z_{\mu }$ fails to adequately approximate $Z_{\text{true}}$ : if $\mu$ is too small, the measurement error will be magnified, whereas an excessively large $\mu$ will lead to an overly smoothed $Z_{\mu }$ , resulting in the loss of local surface details.

3.3. Iterative Compensation

Huang [10] devised a methodology for optimizing the surface reconstruction process by iteratively determining the compensation height value through the gradient residual. The rationale behind this approach lies in the fact that conventional models face difficulties in adequately fitting the high-order items present in the gradient residual data. By incorporating the residual gradient data, compensation for the error induced by these high-order items becomes feasible. Through multiple iterations of compensation, the accuracy of the surface reconstruction can be achieved independently of any underlying model assumptions. The following steps outline the procedure employed in this method:

Step 1: Initialize the improved Southwell model with the initial gradient values $g_0^x$ and $g_0^y$ as presented in this study. Consequently, obtain the initial surface height $z_0\left(x\mathrm{,}y\right)$ and set the current height as $z\left(x\mathrm{,}y\right)=z_0\left(x\mathrm{,}y\right)$ .

Step 2: Calculate the gradients $g^x$ and $g^y$ of the current height $z\left(x\mathrm{,}y\right)$ . Subsequently, derive the residual errors of the x-directional gradient as $g_r^x=g_0^x-g^x$ and the y-directional gradient as $g_r^y=g_0^y-g^y$ .

Step 3: Integrate the obtained gradient residuals to determine the compensation surface height $z_c\left(x\mathrm{,}y\right)$ .

Step 4: Update the current height by incrementing it with $z_c\left(x\mathrm{,}y\right)/n_k$ , where $n_k$ is an empirical parameter, and k represents the iteration number.

Step 5: Repeat Steps 2 to 4 until the height compensation term $z_c\left(x\mathrm{,}y\right)/n_k$ falls below the designated threshold value $z_{thr}$ or the specified number of iterations is reached.

Step 6: Record the resulting height $z\left(x\mathrm{,}y\right)$ as the final surface height.

Huang [10] determined a set of empirical parameters, $n_k=\left[\mathrm{3.0000,4.0909,4.9476,5.6768}\right]$ , through simulations involving surfaces with varying high-order components.

4. Simulation

In the experimental phase, we conducted a comparative analysis between the algorithm proposed in this study and existing algorithms. We denote the Fourier integral method as FC (Frankot Chellappa) [11] , the traditional Southwell model as TFLI (Traditional Finite-difference-based Least-squares Integration) [6] , the method proposed by Li et al. as HFLI (Higher-order Finite-difference-based Least-squares Integration) [7] , and the method proposed by Huang et al. as SLI (Spline-based Least-squares Integration) [9] . Moreover, we introduce the algorithm presented in this paper as ALI (Adams-based Least-square Integration), which does not employ iterative compensation, and label it as Non-iter ALI.

The simulation experiments will be conducted using Matlab R2021b within an experimental environment consisting of an Intel Core i5-7300HQ processor with a main frequency of 2.5 GHz, 20 GB of memory, and a Windows 10 operating system. The simulation experiment follows the subsequent procedure:

1) The surface to be constructed will be defined by Equation (14), encompassing both high-frequency and low frequency information.

2) By calculating the first-order partial derivative of the height data derived from Equation (14), the gradient distribution of the ideal surface in the x-direction and y-direction will be obtained.

3) The gradient data of the ideal surface will serve as input for the FC, TFLI, HFLI, SLI, Non-iter ALI, and ALI algorithms. The resulting height values will be compared against the height data of the ideal surface.

4) The accuracy of the reconstructed surface will be assessed by comparing the RMSE, SSIM, and PV values obtained from each algorithm.

5) The resistance to noise will be evaluated by introducing Gaussian noise of varying magnitudes to the gradient data, and the RMSE value of the height data of the ideal surface will be compared.

6) The computational efficiency will be compared by examining the time required to calculate and execute surface reconstructions of varying sizes.

$z=0.3x\sin \left(x\right)\sin \left(y\right).$ (14)

Among them, the ideal surface to be built selects the surface represented by (14), and its standard figure is shown in Figure 2.

Figure 2. Ideal surface to be built.

By computing the first-order partial derivative of formula (14), we can analyze the distribution of the gradient in both the x and y directions, as illustrated in

Figure 3. Here, we denote the partial derivative with respect to x as $p=\frac{\text{d}z}{\text{d}x}$ and the partial derivative with respect to y as $q=\frac{\text{d}z}{\text{d}y}$ . Notably, it becomes apparent

that significant variations occur in the gradient when there are sharp changes in the morphology.

4.1. Algorithm Precision

Set the size of the reconstruction surface to 128*128, and use the root mean square error RMSE, structural similarity SSIM and peak-to-valley PV to evaluate the reconstruction quality of different algorithms.

RMSE measures the deviation between the predicted value and the true value, expressed as follows:

$\text{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\left(Z_{\text{true}}-Z\right)}^2}$ (15)

Here: $Z_{\text{true}}$ represents the height value of the ideal surface.

SSIM is an index to measure the similarity between two images, which is expressed as follows:

$\text{SSIM}\left(x\mathrm{,}y\right)=\frac{\left(2{\mu }_x{\mu }_y+c_1\right)\left(2{\sigma }_{xy}+c_2\right)}{\left({\mu }_x^2+{\mu }_y^2+c_1\right)\left({\sigma }_x^2+{\sigma }_y^2+c_2\right)}$ (16)

Here: $\mu$ represents the mean value, $\sigma$ represents the standard deviation, $c_1\mathrm{,}{c}_2$ represent the covariance.

PV measures the difference between the maximum value and the minimum value in the surface shape error, expressed as follows:

$\text{PV}=W_{\max }-W_{\min }$ (17)

Here: W represents the error matrix.

The reconstruction surface evaluation results are shown in Table 1.

The comparative analysis depicted in Table 1 unmistakably reveals the algorithm ALI’s superior reconstructive precision, surpassing other algorithms by orders of magnitude. Non-iterative ALI exhibits an approximately 12.2%

Figure 3. Gradient distribution map in x and y direction.

Table 1. Comparison of RMSE, SSIM and PV of each algorithm.

diminution in root mean square error (RMSE) compared to TFLI, whereas ALI showcases a substantial 88.6% reduction in RMSE relative to TFLI. Furthermore, ALI demonstrates notable advancements in terms of structural similarity index (SSIM) and point-to-vertex (PV) evaluation metrics.

Figure 4 illustrates the residual analysis outcomes, underscoring that pronounced surface undulations significantly undermine reconstruction accuracy, especially pertaining to edge shape. In the general Southwell model, the distribution of errors in the surrounding edge region manifests more prominently than in the central area. However, the utilization of the ALI algorithm, incorporating an iterative compensation approach, greatly mitigates this issue. Notably, ALI outperforms other algorithms in 3D reconstruction tasks, particularly for objects with steep topography.

4.2. Noise Immunity

To assess the resilience of each algorithm against noise, Gaussian noise is introduced to the pristine gradient data of size 128*128. Specifically, the mean is established as 0, while the standard deviation is set at 5%, 10%, 15%, and 20%. Figure 5 exhibits a comparative analysis of the anti-noise capabilities of each algorithm.

The observations from Figure 5 reveal a noteworthy trend: as the standard deviation of Gaussian noise escalates, the root mean square error (RMSE) values exhibit an upward trajectory for all examined techniques. However, the RMSE of the proposed ALI algorithm in this study remains notably lower than alternative methods. Even in the absence of iterative compensation using Non-iter ALI, the

Figure 4. Comparison chart of residual errors on the surface to be built.

Figure 5. Comparison of noise resistance of various algorithms.

incorporation of Tikhonov Regularization ensures a diminished error when compared to other Southwell model approaches. Furthermore, with increasing noise, the reconstructed surface generated by the ALI method gradually adopts a rougher texture. Nevertheless, it still manages to capture a highly accurate representation of the surface shape. This exemplifies the robust noise resilience of the algorithm outlined in this research endeavor.

4.3. Time Consuming

Apart from precision and resistance to noise, computational time is also a crucial factor to be considered. Notably, the time-consuming nature of calculations varies as the reconstructed surface size fluctuates across different algorithms. To provide a comprehensive assessment of the computational time requirements for the algorithm proposed in this study, we conducted evaluations with surface sizes set at 128128, 256256, and 512*512, respectively. The outcomes of these evaluations are illustrated in Figure 6.

Figure 6. Calculation time chart of each algorithm.

The provided figure elucidates pertinent insights regarding the reconstruction speed of different methods. Remarkably, the FC method exhibits the swiftest reconstruction pace, and its increment in time consumption remains minimal as the reconstructed surface expands. Among the three methods-TFLI, HFLI, and SLI-they demonstrate comparable time consumption for reconstruction across various surface sizes. Notably, the proposed Non-iter ALI method showcases lower time consumption than these three methods. While the AIL method excels in accuracy and noise resilience, its reliance on the iterative compensation approach leads to the lengthiest reconstruction time, approximately five times that of Non-iter ALI. Consequently, the selection of an appropriate algorithm should be guided by the characteristics of the reconstructed surface and the desired reconstruction requirements.

5. Conclusion and Suggestion

This paper employs the fourth-order Adams method to enhance the connection between the measured gradient and unknown height in the conventional Southwell model. Simultaneously, the algorithm incorporates Tikhonov regularization and iterative compensation techniques to bolster anti-noise capabilities and reconstruction accuracy. By contrasting the performance of the Non-iter ALI, ALI, and FC methods proposed herein with other Southwell model approaches, the experimental results unveil the following findings: 1) Under ideal conditions, ALI surpasses other methods by order of magnitude in reconstruction accuracy. 2) Both Non-iter ALI and ALI exhibit superior anti-noise capabilities compared to alternative algorithms in varying noise environments. 3) Non-iter ALI boasts the shortest reconstruction time when compared to TFLI, HFLI, and SLI methods, while ALI experiences the lengthiest reconstruction time, approximately five times that of Non-iter ALI.

Fund Project

This work was supported by the Jishou University Graduate Research Innovation Project (JGY2023072).

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Li, Y., Sixou, B. and Peyrin, F. (2021) A Review of the Deep Learning Methods for Medical Images Super Resolution Problems. IRBM, 42, 120-133. https://doi.org/10.1016/j.irbm.2020.08.004
[2]	Mochi, I. and Goldberg, K.A. (2015) Modal Wavefront Reconstruction from Its Gradient. Applied Optics, 54. 3780-3785. https://doi.org/10.1364/AO.54.003780
[3]	Ghosh, A., Fyffe, G., Tunwattanapong, B., Busch, J., Yu, X. and Debevec, P. (2011) Multiview Face Capture Using Polarized Spherical Gradient Illumination. ACM Transactions on Graphics (TOG), 30, 1-10. https://doi.org/10.1145/2070781.2024163
[4]	Hudgin, R.H. (1977) Optimal Wave-Front Estimation. Journal of the Optical Society of America, 67, 378-382. https://doi.org/10.1364/JOSA.67.000378
[5]	Fried, D.L. (1977) Least-Square Fitting A Wave-Front Distortion Estimate to An Array of Phase-Difference Measurements. Journal of the Optical Society of America, 67, 370-375. https://doi.org/10.1364/JOSA.67.000370
[6]	Southwell, W.H. (1980) Wave-Front Estimation from Wave-Front Slope Measurements. Journal of the Optical Society of America, 70, 998-1006. https://doi.org/10.1364/JOSA.70.000998
[7]	Li, G., Li, Y., Liu, K., Ma, X. and Wang, H. (2013) Improving Wavefront Reconstruction Accuracy by Using Integration Equations with Higher-Order Truncation Errors in the Southwell Geometry. Journal of the Optical Society of America A, 30, 1448-1459. https://doi.org/10.1364/JOSAA.30.001448
[8]	Ren, H., Gao, F. and Jiang, X. (2016) Least-Squares Method for Data Reconstruction from Gradient Data in Deflectometry. Applied Optics, 55, 6052-6059. https://doi.org/10.1364/AO.55.006052
[9]	Huang, L., Xue, J., Gao, B., Zuo, C. and Idir, M. (2017) Spline Based Least Squares Integration for Two-Dimensional Shape or Wavefront Reconstruction. Optics and Lasers in Engineering, 91, 221-226. https://doi.org/10.1016/j.optlaseng.2016.12.004
[10]	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
[11]	Frankot, R. T and Chellappa, R. (1988) A Method for Enforcing Integrability in Shape from Shading Algorithms. IEEE Transactions on Pattern Analysis and Machine Intelligence, 10, 439-451. https://doi.org/10.1109/34.3909