Acquisition, creation and processing of 3D digital objects is an important topic in various fields, e.g., medical imaging , computer graphics [11, 8], industry , and preservation of cultural heritage . A fundamental step is to reconstruct a surface from a 3D scanned point cloud data , denoted as for or , such as in Figure 1.
We focus on reconstructing an -dimensional manifold , a curve in or a surface in , from the point cloud
. We assume only the point locations are given, and no other geometrical information such as normal vectors at each point are known. We explore fast algorithms for minimizing the following energy proposed in:
where is the distance from any point to , is a positive integer, and is the surface area element. This energy is the -norm of the distance function restricted in . In , the authors used the fast sweeping scheme to solve Eikonal equation, and the Euler-Lagrange equations are solved by gradient descent algorithm.
There are various related works on surface reconstruction from point cloud data: a convection model proposed in , a data-driven logarithmic prior for noisy data in , using surface tension to enrich the Euler-Lagrange equations in 
, and using principal component analysis to reconstruct curves embedded in sub-manifolds in. A semi-implicit scheme is introduced in  to simulate the curvature and surface diffusion motion of the interface. In , the authors defined the surface via a collection of anisotropic Gaussians centered at each entry of the input point cloud, and used TVG-L1 model  for minimization. A similar strategy addresses an gradient regularization model proposed in . Some models incorporate additional information. In , the authors proposed a novel variational model, consisting of the distance, the normal, and the smoothness term. Euler’s Elastica model is incorporated for surface reconstruction in  where graph cuts algorithm is used. The model in  extends the active contours segmentation model to 3D and implicitly allows controlling the curvature of the level set function.
In this paper, we explore fast algorithms to minimize the weighted minimum surface energy (1) for and 2. We propose a Semi-Implicit Method (SIM) to relax the time-step constraint for , and an Augmented Lagrangian Method (ALM) based on the alternating direction method of multipliers (ADMM) approach for . These algorithms minimize the weighted minimal surface energy (1) with high accuracy and superior efficiency. We analyze the behavior of ALM in terms of parameter choices and explore its connection to SIM. Various numerical experiments are presented to discuss the effects of the algorithms.
2 Proposed Algorithms
Let ( or ) denote a bounded domain containing the given -dimensional point cloud data, . Using the level-set formulation for , the -weighted minimum surface energy (1) can be rewritten as:
Here is the Dirac delta function which takes when , and elsewhere. Compared to (1), this integral is defined on , which makes the computation flexible and free from explicitly tracking . We use for SIM introduced in Section 2.1, and for ALM in Section 2.2. In general, is a natural choice, since it provides better stability and efficiency for a semi-implicit type PDE-based method. For ALM, we explore to take advantage of an aspect of fast algorithm in ADMM setting such as shrinkage, similarly to the case in . Visually, the numerical results of surface reconstruction are similar for or (see Section 3).
2.1 Semi-Implicit Method (SIM) to minimize
We introduce a gradient-flow-based semi-implicit method to Minimize
Following , the first variation of with respect to is characterized as a functional:
for any test function from the Sobolev space . Minimizing (3) is equivalent to finding the critical point such that This is associated with solving the following initial value problem:
where is an initial guess for the unknown , and The steady state solution of (4) gives a minimizer of .
Here the delta function is realized as the derivative of the one dimensional Heaviside function . We adopt the smooth approximation of as in :
with as the smoothness parameter. Then is approximated by its smoothed version expressed as
Employing a semi-implicit scheme, we solve from (6) by iteratively updating using via the following equation:
is the time-step. This equation can be efficiently solved by the Fast Fourier Transform (FFT). Denoting the discrete Fourier transform byand its inverse by , we have:
Accordingly, the discrete Fourier transform of is
As for the stopping criterion, we exploit the mean relative change of the weighted minimum surface energy (1). At the iteration, the algorithm terminates if
Here the quantity represents the average of the energy values computed from the to the iteration for some , . We fix , and set for SIM. We summarize the main steps of SIM in Algorithm 1.
2.2 Augmented Lagrangian Method (ALM) to Minimize
In this section, we present an augmented Lagrangian-based method to minimize the weighted minimum surface energy (2) for , i.e.,
For the non-differentiable term in (10), we utilize the variable-splitting and introduce an auxiliary variable . We rephrase the minimization of as a constrained optimization problem:
where is a scalar penalty parameter and represents the Lagrangian multiplier. Minimizing (12) amounts to considering the following saddle-point problem:
Given , , and , for , the iteration of an ADMM-type algorithm for (13) consists of solving a series of sub-problems:
Each sub-problem can be solved efficiently. First, we find the minimizer of the sub-problem (14) by solving its Euler-Lagrange equation:
Here is the Laplacian operator, and we solve this via FFT, similar to (8) for SIM. Thus, the sub-problem is solved via:
Second, the sub-problem (15) is equivalent to a weighted Total Variation (TV) minimization, whose solution admits a closed-form expression using the shrinkage operator . Explicitly, the updated is computed via:
2.3 Connection between SIM and ALM algorithms
Note that both SIM and ALM involve solving elliptic PDEs of the form:
for some constants , and a function defined on . For SIM, it is equation (7):
and for ALM, it is equation (18):
We remark interesting connections between SIM and ALM. First, both methods have stabilizing terms but in different positions on the left side of (21). For SIM, it is , while for ALM, it is . Second, relating the coefficients of , in SIM gives insight to the effect of in ALM. In general, a large slows down the convergence of ALM, while a small accelerates it (as the effect of on SIM). Figure 2 shows convergence behaviors of ALM for different , using the five-fold circle point cloud in Figure 1 (a). It displays the CPU time (in seconds) for , , and varying from to . Note that as increases, the time required to reach the convergence increases almost quadratically at first, then stays around the same level. Third, the correspondence between in SIM, and in ALM allows another interpretation of the parameter in ALM. In SIM, a large smears the solution and avoids discontinuities or sharp corners, and for ALM, large also allows to pass through fine details. Figure 7 in Section 3 presents more details, where we experiment with different and values for the five-fold circle point cloud shown in Figure 1 (a).
3 Numerical Implementations, Experiments and Effects of Parameters
In this section, we describe the implementation details and present numerical experiments. For both SIM and ALM, we vary from to . For SIM, we use . When is in 2D, we set , and for 3D. For ALM, ranges from to , and from to .
The code is written in Matlab and executed without additional machine support, e.g. parallelization or GPU-enhanced computations. All the experiments are performed on Intel® Core™4-Core 1.8GHz (4.0GHz with Turbo) machine, with 16 GB/RAM and Intel® UHD Graphics 620 graphic card under Windows OS. The contours and isosurfaces are displayed using Matlab visualization engine. No post-processing, e.g., smoothing or sharpening, is applied.
3.1 Implementation Details
We illustrate the details for planar point clouds, i.e., , and the extension to is straightforward. Let the computational domain , , be discretized by a Cartesian grid with . For any function (or a vector field ) defined on , we use or to denote . We use the usual backward and forward finite difference schemes:
The gradient, divergence and the Laplacian operators are approximated as follows:
The distance function is computed once at the beginning and no update is needed. It satisfies an Eikonal equation:
and discretizing (22) via the Lax-Friedrich scheme leads to an updating formula:
Keeping to be a signed distance function during the iteration improves the stability of level-set-based algorithms. We reinitialize at the iteration by solving the following PDE:
Here subscript represents artificial time and partial derivative, and is the sign function. In practice, being a signed distance function near -level-set is important; thus, it is sufficient to solve (24) for a few steps. We fix 10 steps of reinitialization throughout this paper.
3.2 Numerical Experiments of 2D and 3D Point Clouds
The first experiment, Figure 3, is reconstruction of planar curves from 2D point clouds confined within a square . We generate the data using four different shapes: a triangle, an ellipse, a square whose corners are missing, and a five-fold circle. For these cases, we use a centered circle with radius as the initial guess, shown in Figure 3 (a). Figure 3 (b) and (c) display the given , as well as the curves identified by SIM and ALM with , respectively. Both methods produce comparably accurate results. In the triangle example, corners get as close as the approximated delta function (with parameter ) allows for both methods. The ellipse and square results fit very closely to the given point clouds. For the five-fold-circle, there is a slight difference in how the curve fits the edges, yet the results are very compatible.
Table 1 shows the CPU time (in seconds) for SIM, ALM using , and , and the explicit method in  with on the same data sets. With proper choices of , ALM outperforms the other methods in terms of computational efficiency. SIM is stable without any dependency on the choice of parameters, and its run-times are comparable to the best performances of ALM in most cases. Both methods are faster than the explicit method in all the examples, and computational efficiency is similar between ALM and SIM.
The second set of experiments reconstruct surface from the point cloud in 3D: a jar in Figure 1 (b) and a torus in Figure 1 (c) within . In Figure 4, we show the reconstructed surfaces using SIM and ALM. A portion of the given point cloud is superposed for validation in each case. Both methods successfully capture the overall shapes and non-convex features of the jar, as well as the torus. There are only slight differences in the reconstruction between using SIM with and using ALM with .
Table 2 shows the efficiency of SIM and ALM compared to the explicit method in  for the experiments in Figure 4. Thanks to the semi-implicit scheme, can be large and we used in SIM; in the explicit method, we are forced to use much smaller time step to maintain the stability. The improvement of run-time in ALM is carefully controlled by the parameters , and . We choose , and for both cases. Both SIM and ALM efficiently provide accurate reconstruction.
The third set of examples show the effect of the distance function . Notice that the weighted minimal surface energy (1) is mainly driven by the distance function , that is, the given point cloud determines the landscape of , which affects the behavior of the level-set during the evolution. Figure 5 shows the evolution using ALM, applied to different subsets of point clouds sampled from the same bunny face shape. The densities of the point cloud vary for the three different regions: the face with points, the head with points, and each ear with points. Figure 5 (a) shows the given point cloud for , with the -level-set of at iteration, (b) for , at iteration, and(c) for , at iteration. These three curves eventually degenerate to a point. (d) for and shows the converged solution. In (a)–(c), denser parts of the point cloud attract the curve with stronger forces, the sparser parts of the point cloud fail to lock the curve. Then, the energy model (2) drives curves to have short lengths, i.e., the level set tends to shrink. In (d), with a more balanced distribution of points, the curve converges to correct shape.
The fourth set of examples demonstrate the robustness of ALM and SIM against noise. Figure 6
shows the reconstructed curves from clean and noisy data: (a)-(c) are results of ALM, and (d)-(f) are results of SIM. (a) and (d) in the first column show results obtained from the clean data, which has 200 points sampled from a three-fold circle. Gaussian noise with standard deviationis added to both and coordinates to generate noisy point cloud in the second column, (b) and (e). To show the differences, the third column superposes both results reconstructed from clean and noisy point clouds. Both ALM and SIM provide compatible results. For the noisy data, although the reconstructed curves show some oscillation, they are very close to the solutions using the clean data, respectively.
3.3 Choice of Parameters for ALM and its Effects
The proposed ALM has one parameter , and the model (2) uses the delta function, where the smoothness parameter is added to stabilize the computation. Both parameters have straightforward effects on the level-set evolution from (17). For example, consider a set of points within a thin-band around the 0-level-set of , denoted by . By the continuity of , there exist and such that and ; these values are the minimum and maximum of the function , respectively. At these points, (17) takes the following forms:
The first terms in the right hand side of (25) show that with a smaller , there are less number of points in , but influence from is stronger. With a larger , affects more number of points in , but the influence becomes weaker. Varying also modifies the effect of , while the size of is not changed.
We also find that interacts with and effectively modifies the shape of the level-set. Figure 7 shows the results for ALM using different combinations of and , on the five-fold circle point cloud in Figure 1 (a). For a fixed , increasing makes the approximated delta function smoother; consequently, narrow and elongated shapes are omitted, and the reconstructed curve becomes more convex. For a fixed , larger loses more details, as discussed in Subsection 2.3. The speed of convergence varies for different combinations of and . When the choices are reasonable, the algorithm converges fast within 2 seconds. When both and are large, results are not as good, and the convergences are slow.
Another observation comes from (20). For any point and , if the following value:
is positive, then , and has no direct effect on (18) at in the next iteration. Regarding as a quadratic polynomial in terms of parameterized by , the sign of depends on the sign of and the sign of its discriminant computed via:
The sign of is related to the position of relative to the 0-level-set. The sign of is determined by comparing the length of a vector difference with the quantity . By the projection theorem, is bounded below by , i.e., the squared residual of orthogonal projection of onto ; therefore, we can decide the sign of using via the following cases:
When , for any , .
if or , then ;
if , then .
When , concaves upwards and for any . If , is positive for all and has no effect on level set evolution. If , is positive for outside the interval bounded by two roots of , i.e.,
When , concaves downwards, and for any . In this case, is never positive: either , i.e., no roots, or but both roots are negative.
Notice that the bounds, and , are closely related to the ratio , which contributes to the adaptive behavior of ALM. For example, for a point where , when is close to but , and becomes extremely large; thus, for a moderate value of , has strong influence on the evolution of the level-set near and swiftly moves the curve towards the point cloud. For a point which is close to both and , the level-set evolution becomes more stringent about the minimization of the energy (10).
Figure 8 illustrates this effect, for the five-fold circle point cloud in Figure 1 (a) with and . Figure 8 shows (a) , (b) , (c) , and (d) the region where effects the level set evolution. The figures are for iterations and (converged). The region inside always experiences the influence of , as described above. Figure 8 (a) shows that the region outside is mostly blue indicating ; hence, for almost every point outside the 0-level-set, as long as , the landscape of has strong effects on the evolution. In (b) and (c), observe that high values of only concentrate near the 0-level-set while remains relatively small in the whole domain; thus, the influence of is strong near . (d) displays the white regions where explicitly guides the level-set evolution and the black regions where has no direct effect. These results show that, although ALM evolves the level-set globally, it ignores the effects of when evolving the regions far away from the level-sets; and it utilizes to refine the local structures for the regions of the level-sets close to .
We propose two fast algorithms, SIM and ALM, to reconstruct -dimensional manifold from unstructured point clouds by minimizing the weighted minimum surface energy (2). SIM improves the computational efficiency by relaxing the constraint on the time-step using a semi-implicit scheme. ALM follows an augmented Lagrangian approach and solves the problem by an ADMM-type algorithm. Numerical experiments show that the proposed algorithms are superior at the computational speed, and both of them produce accurate results. Theoretically, we demonstrate the delicate interaction among parameters involved in ALM, and show the connections between SIM and ALM. This explains the behaviors of ALM from the perspective of SIM.
-  M. Alexa, J. Behr, D. Cohen-Or, S. Fleishman, D. Levin, and C. T. Silva. Point set surfaces. In Proceedings of the Conference on Visualization’01, pages 21–28. IEEE Computer Society, 2001.
-  M. Alexa, J. Behr, D. Cohen-Or, S. Fleishman, D. Levin, and C. T. Silva. Computing and rendering point set surfaces. IEEE Transactions on Visualization and Computer Graphics, 9(1):3–15, January 2003.
-  E. Bae, X.-C. Tai, and W. Zhu. Augmented Lagrangian method for an Euler’s elastica based segmentation model that promotes convex contours. Inverse Problems & Imaging, 11(1):1–23, 2017.
-  Z. Bi and L. Wang. Advances in 3D data acquisition and processing for industrial applications. Robotics and Computer-Integrated Manufacturing, 26(5):403–413, 2010.
-  M. Bolitho, M. Kazhdan, R. Burns, and H. Hoppe. Parallel poisson surface reconstruction. In International symposium on visual computing, pages 678–689. Springer, 2009.
-  R. Bracewell and R. Bracewell. The Fourier Transform and Its Applications. Electrical Engineering Series. McGraw Hill, 2000.
-  X. Bresson, S. Esedoḡlu, P. Vandergheynst, J.-P. Thiran, and S. Osher. Fast global minimization of the active contour/snake model. Journal of Mathematical Imaging and vision, 28(2):151–167, 2007.
-  F. Calakli and G. Taubin. SSD: Smooth signed distance surface reconstruction. In Computer Graphics Forum, volume 30, pages 1993–2002. Wiley Online Library, 2011.
-  J. C. Carr, R. K. Beatson, J. B. Cherrie, T. J. Mitchell, W. R. Fright, B. C. McCallum, and T. R. Evans. Reconstruction and representation of 3D objects with radial basis functions. In Proceedings of the 28th Annual Conference on Computer Graphics and Interactive Techniques, pages 67–76. ACM, 2001.
-  J. C. Carr, R. K. Beatson, B. C. McCallum, W. R. Fright, T. J. McLennan, and T. J. Mitchell. Smooth surface reconstruction from noisy range data. In Proceedings of the 1st International Conference on Computer Graphics and Interactive Techniques in Australasia and South East Asia, pages 119–ff. ACM, 2003.
G. Casciola, D. Lazzaro, L. B. Montefusco, and S. Morigi.
Shape preserving surface reconstruction using locally anisotropic radial basis function interpolants.Computers & Mathematics with Applications, 51(8):1185–1198, 2006.
-  T. F. Chan and L. A. Vese. Active contours without edges. IEEE Transactions on Image Processing, 10(2):266–277, 2001.
H. Q. Dinh, G. Turk, and G. Slabaugh.
Reconstructing surfaces using anisotropic basis functions.
Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001, volume 2, pages 606–613. IEEE, 2001.
-  V. Estellers, M. Scott, K. Tew, and S. Soatto. Robust poisson surface reconstruction. In International Conference on Scale Space and Variational Methods in Computer Vision, pages 525–537. Springer, 2015.
-  V. Estellers, D. Zosso, R. Lai, S. Osher, J.-P. Thiran, and X. Bresson. Efficient algorithm for level set method preserving distance function. IEEE Transactions on Image Processing, 21(12):4722–4734, 2012.
-  L. Gomes, O. R. P. Bellon, and L. Silva. 3D reconstruction methods for digital preservation of cultural heritage: A survey. Pattern Recognition Letters, 50:3–14, 2014.
-  J. Haličková and K. Mikula. Level set method for surface reconstruction and its application in surveying. Journal of Surveying Engineering, 142(3):04016007, 2016.
-  H. Huang, D. Li, H. Zhang, U. Ascher, and D. Cohen-Or. Consolidation of unorganized point clouds for surface reconstruction. ACM Transactions on Graphics (TOG), 28(5):176, 2009.
-  C. Y. Kao, S. Osher, and J. Qian. Lax–Friedrichs sweeping scheme for static Hamilton–Jacobi equations. Journal of Computational Physics, 196(1):367–391, 2004.
-  M. Kazhdan, M. Bolitho, and H. Hoppe. Poisson surface reconstruction. In Proceedings of the 4th Eurographics Symposium on Geometry Processing, volume 7, pages 61–70, 2006.
-  M. Kazhdan and H. Hoppe. Screened poisson surface reconstruction. ACM Transactions on Graphics (TOG), 32(3):29, 2013.
-  D. Khan, M. A. Shirazi, and M. Y. Kim. Single shot laser speckle based 3D acquisition system for medical applications. Optics and Lasers in Engineering, 105:43–53, 2018.
-  H. Li, Y. Li, R. Yu, J. Sun, and J. Kim. Surface reconstruction from unorganized points with gradient minimization. Computer Vision and Image Understanding, 169:108–118, 2018.
-  X. Li, W. Wan, X. Cheng, and B. Cui. An improved Poisson surface reconstruction algorithm. In 2010 International Conference on Audio, Language and Image Processing, pages 1134–1138. IEEE, 2010.
-  J. Liang, F. Park, and H.-K. Zhao. Robust and efficient implicit surface reconstruction for point clouds based on convexified image segmentation. Journal of Scientific Computing, 54(2-3):577–602, 2013.
-  H. Liu, X. Wang, and W. Qiang. Implicit surface reconstruction from 3D scattered points based on variational level set method. In 2008 2nd International Symposium on Systems and Control in Aerospace and Astronautics, pages 1–5. IEEE, 2008.
-  H. Liu, Z. Yao, S. Leung, and T. F. Chan. A level set based variational principal flow method for nonparametric dimension reduction on Riemannian manifolds. SIAM Journal on Scientific Computing, 39(4):A1616–A1646, 2017.
-  S. Osher and R. P. Fedkiw. Level set methods: an overview and some recent results. Journal of Computational Physics, 169(2):463–502, 2001.
-  S. Osher and J. A. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. Journal of Computational Physics, 79(1):12–49, 1988.
-  A. C. Öztireli, G. Guennebaud, and M. Gross. Feature preserving point set surfaces based on non-linear kernel regression. In Computer Graphics Forum, volume 28, pages 493–501. Wiley Online Library, 2009.
-  J. A. Sethian. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid mechanics, Computer vision, and Materials Science, volume 3. Cambridge university press, 1999.
-  J. Shi, M. Wan, X.-C. Tai, and D. Wang. Curvature minimization for surface reconstruction with features. In International Conference on Scale Space and Variational Methods in Computer Vision, pages 495–507. Springer, 2011.
-  Y. Shi and W. C. Karl. Shape reconstruction from unorganized points with a data-driven level set method. In 2004 IEEE International Conference on Acoustics, Speech, and Signal Processing, volume 3, pages iii–13. IEEE, 2004.
-  P. Smereka. Semi-implicit level set methods for curvature and surface diffusion motion. Journal of Scientific Computing, 19(1):439–456, 2003.
-  X.-C. Tai, J. Hahn, and G. J. Chung. A fast algorithm for Euler’s elastica model using augmented Lagrangian method. SIAM Journal on Imaging Sciences, 4(1):313–344, 2011.
-  A. Tsai, A. Yezzi Jr, W. Wells, C. Tempany, D. Tucker, A. Fan, W. E. Grimson, and A. Willsky. A shape-based approach to the segmentation of medical imagery using level sets. IEEE Transactions on Medical Imaging, 22(2):137, 2003.
-  H.-K. Zhao, S. Osher, and R. Fedkiw. Fast surface reconstruction using the level set method. In Proceedings IEEE Workshop on Variational and Level Set Methods in Computer Vision, pages 194–201. IEEE, 2001.
-  H.-K. Zhao, S. Osher, B. Merriman, and M. Kang. Implicit, nonparametric shape reconstruction from unorganized points using a variational level set method. Computer Vision and Image Understanding, 80(3):295–319, 2000.