1 Introduction
Acquisition, creation and processing of 3D digital objects is an important topic in various fields, e.g., medical imaging [22], computer graphics [11, 8], industry [4], and preservation of cultural heritage [16]. A fundamental step is to reconstruct a surface from a 3D scanned point cloud data [1], denoted as for or , such as in Figure 1.
(a)  (b)  (c) 

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
[38]:(1) 
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 [38], the authors used the fast sweeping scheme to solve Eikonal equation, and the EulerLagrange equations are solved by gradient descent algorithm.
Among many ways to represent the underlying surface, e.g., (moving) least square projection [2, 30]
[13, 9, 10, 11, 18], poisson reconstruction [20, 5, 24, 21, 14], we use the level set method as in [15, 37, 38]. The level set formulation allows topological changes and selfintersection during the evolution [29, 28], and gain popularity in many applications [12, 36, 31] just to mention a few. We represent the surface as a zero level set of :There are various related works on surface reconstruction from point cloud data: a convection model proposed in [37], a datadriven logarithmic prior for noisy data in [33], using surface tension to enrich the EulerLagrange equations in [17]
, and using principal component analysis to reconstruct curves embedded in submanifolds in
[27]. A semiimplicit scheme is introduced in [34] to simulate the curvature and surface diffusion motion of the interface. In [25], the authors defined the surface via a collection of anisotropic Gaussians centered at each entry of the input point cloud, and used TVGL1 model [7] for minimization. A similar strategy addresses an gradient regularization model proposed in [23]. Some models incorporate additional information. In [26], 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 [32] where graph cuts algorithm is used. The model in [15] 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 SemiImplicit Method (SIM) to relax the timestep 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 levelset formulation for , the weighted minimum surface energy (1) can be rewritten as:
(2) 
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 semiimplicit type PDEbased 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 [3]. Visually, the numerical results of surface reconstruction are similar for or (see Section 3).
2.1 SemiImplicit Method (SIM) to minimize
We introduce a gradientflowbased semiimplicit method to Minimize
(3) 
Following [38], 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:
(4) 
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 [6]:
(5) 
with as the smoothness parameter. Then is approximated by its smoothed version expressed as
We add a stabilizing diffusive term for on both sides of the PDE in (4) to consolidate the computation, similarly to [34]:
(6) 
Employing a semiimplicit scheme, we solve from (6) by iteratively updating using via the following equation:
(7) 
where
is the timestep. This equation can be efficiently solved by the Fast Fourier Transform (FFT). Denoting the discrete Fourier transform by
and its inverse by , we have:Accordingly, the discrete Fourier transform of is
Here the coefficient in front of represents the diagonalized discrete Laplacian operator in the frequency domain. Let be the right side of (7), then the solution of (7) is computed via
(8) 
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
(9) 
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 Lagrangianbased method to minimize the weighted minimum surface energy (2) for , i.e.,
(10) 
For the nondifferentiable term in (10), we utilize the variablesplitting and introduce an auxiliary variable . We rephrase the minimization of as a constrained optimization problem:
(11) 
here we replace by its smooth approximation as in (5). To solve problem (11), we formulate the augmented Lagrangian function:
(12) 
where is a scalar penalty parameter and represents the Lagrangian multiplier. Minimizing (12) amounts to considering the following saddlepoint problem:
(13) 
Given , , and , for , the iteration of an ADMMtype algorithm for (13) consists of solving a series of subproblems:
(14)  
(15)  
(16) 
Each subproblem can be solved efficiently. First, we find the minimizer of the subproblem (14) by solving its EulerLagrange equation:
(17) 
Following [3], we introduce a frozencoefficient term , for , on both sides of (17) to stabilize the computation; thus, (14) is solved using the following equation:
(18) 
Here is the Laplacian operator, and we solve this via FFT, similar to (8) for SIM. Thus, the subproblem is solved via:
(19) 
2.3 Connection between SIM and ALM algorithms
Note that both SIM and ALM involve solving elliptic PDEs of the form:
(21) 
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 fivefold 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 fivefold 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 GPUenhanced computations. All the experiments are performed on Intel^{®} Core™4Core 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 postprocessing, 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:
(22) 
and discretizing (22) via the LaxFriedrich scheme leads to an updating formula:
(23) 
We solve (23) using the fast sweeping method [19] with complexity for grid points.
Keeping to be a signed distance function during the iteration improves the stability of levelsetbased algorithms. We reinitialize at the iteration by solving the following PDE:
(24) 
Here subscript represents artificial time and partial derivative, and is the sign function. In practice, being a signed distance function near levelset 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 fivefold 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 fivefoldcircle, 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 [38] 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 runtimes 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.
Object  ALM()  ALM(  ALM()  ALM()  SIM  [38] 
Triangle  
Ellipse  
Square  1.20  
Fivefold circle 
(a)  (b) 
(c)  (d) 
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 nonconvex 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 [38] for the experiments in Figure 4. Thanks to the semiimplicit 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 runtime in ALM is carefully controlled by the parameters , and . We choose , and for both cases. Both SIM and ALM efficiently provide accurate reconstruction.
Object  ALM  SIM  [38] 

Jar  
Torus 
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 levelset 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 levelset 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 threefold circle. Gaussian noise with standard deviation
is 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.(a)  (b)  (c) 
(d)  (e)  (f) 
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 levelset evolution from (17). For example, consider a set of points within a thinband around the 0levelset 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:
(25) 
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 levelset. Figure 7 shows the results for ALM using different combinations of and , on the fivefold 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 0levelset. 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 , .

When :

if or , then ;

if , then .
Here,

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 levelset near and swiftly moves the curve towards the point cloud. For a point which is close to both and , the levelset evolution becomes more stringent about the minimization of the energy (10).
(a)  (b) 
(c)  (d) 
Figure 8 illustrates this effect, for the fivefold 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 0levelset, 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 0levelset while remains relatively small in the whole domain; thus, the influence of is strong near . (d) displays the white regions where explicitly guides the levelset evolution and the black regions where has no direct effect. These results show that, although ALM evolves the levelset globally, it ignores the effects of when evolving the regions far away from the levelsets; and it utilizes to refine the local structures for the regions of the levelsets close to .
4 Conclusion
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 timestep using a semiimplicit scheme. ALM follows an augmented Lagrangian approach and solves the problem by an ADMMtype 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.
References
 [1] M. Alexa, J. Behr, D. CohenOr, 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.
 [2] M. Alexa, J. Behr, D. CohenOr, 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.
 [3] 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.
 [4] Z. Bi and L. Wang. Advances in 3D data acquisition and processing for industrial applications. Robotics and ComputerIntegrated Manufacturing, 26(5):403–413, 2010.
 [5] M. Bolitho, M. Kazhdan, R. Burns, and H. Hoppe. Parallel poisson surface reconstruction. In International symposium on visual computing, pages 678–689. Springer, 2009.
 [6] R. Bracewell and R. Bracewell. The Fourier Transform and Its Applications. Electrical Engineering Series. McGraw Hill, 2000.
 [7] 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.
 [8] F. Calakli and G. Taubin. SSD: Smooth signed distance surface reconstruction. In Computer Graphics Forum, volume 30, pages 1993–2002. Wiley Online Library, 2011.
 [9] 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.
 [10] 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.

[11]
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.  [12] T. F. Chan and L. A. Vese. Active contours without edges. IEEE Transactions on Image Processing, 10(2):266–277, 2001.

[13]
H. Q. Dinh, G. Turk, and G. Slabaugh.
Reconstructing surfaces using anisotropic basis functions.
In
Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001
, volume 2, pages 606–613. IEEE, 2001.  [14] 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.
 [15] 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.
 [16] 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.
 [17] 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.
 [18] H. Huang, D. Li, H. Zhang, U. Ascher, and D. CohenOr. Consolidation of unorganized point clouds for surface reconstruction. ACM Transactions on Graphics (TOG), 28(5):176, 2009.
 [19] 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.
 [20] 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.
 [21] M. Kazhdan and H. Hoppe. Screened poisson surface reconstruction. ACM Transactions on Graphics (TOG), 32(3):29, 2013.
 [22] 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.
 [23] 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.
 [24] 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.
 [25] 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(23):577–602, 2013.
 [26] 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.
 [27] 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.
 [28] S. Osher and R. P. Fedkiw. Level set methods: an overview and some recent results. Journal of Computational Physics, 169(2):463–502, 2001.
 [29] S. Osher and J. A. Sethian. Fronts propagating with curvaturedependent speed: algorithms based on Hamilton–Jacobi formulations. Journal of Computational Physics, 79(1):12–49, 1988.
 [30] A. C. Öztireli, G. Guennebaud, and M. Gross. Feature preserving point set surfaces based on nonlinear kernel regression. In Computer Graphics Forum, volume 28, pages 493–501. Wiley Online Library, 2009.
 [31] 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.
 [32] 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.
 [33] Y. Shi and W. C. Karl. Shape reconstruction from unorganized points with a datadriven level set method. In 2004 IEEE International Conference on Acoustics, Speech, and Signal Processing, volume 3, pages iii–13. IEEE, 2004.
 [34] P. Smereka. Semiimplicit level set methods for curvature and surface diffusion motion. Journal of Scientific Computing, 19(1):439–456, 2003.
 [35] 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.
 [36] A. Tsai, A. Yezzi Jr, W. Wells, C. Tempany, D. Tucker, A. Fan, W. E. Grimson, and A. Willsky. A shapebased approach to the segmentation of medical imagery using level sets. IEEE Transactions on Medical Imaging, 22(2):137, 2003.
 [37] 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.
 [38] 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.
Comments
There are no comments yet.