Multi-Reconstruction from Points Cloud by Using a Modified Vector-Valued Allen–Cahn Equation

The Poisson surface reconstruction algorithm has become a very popular tool of reconstruction from point clouds. If we reconstruct each region separately in the process of multi-reconstruction, then the reconstructed objects may overlap with each other. In order to reconstruct multicomponent surfaces without self-intersections, we propose an efficient multi-reconstruction algorithm based on a modified vector-valued Allen–Cahn equation. The proposed algorithm produces smooth surfaces and closely preserves the original data without self-intersect. Based on operator splitting techniques, the numerical scheme is divided into one linear equation and two nonlinear equations. The linear equation is discretized using an implicit method, and the resulting discrete system of equation is solved by a fast Fourier transform. The two nonlinear equations are solved analytically due to the availability of a closed-form solution. The numerical scheme has merit in that it can be straightforwardly applied to a graphics processing unit, allowing for accelerated implementation that performs much faster than central processing unit alternatives. Various experimental, numerical results demonstrate the effectiveness and robustness of the proposed method.


Introduction
Multimaterial reconstruction from multi-labeled point clouds is of significant importance in many disciplines, such as mechanical engineering [1,2]. Such point clouds are sampled by a multi-camera vision system and provide the local connectivity of objects. For instance, in the reversed engineering problem, multicomponent surfaces/volumes of the object should be reconstructed without self-intersections so that a multimaterial 3D printer can process these models. The challenges in multi-reconstruction mainly lie in the following two aspects: how to handle the contacting regions and how to maintain the original topology. Up until now, many reconstruction methods have been developed to address modeling from point clouds. In addition, modeling multi-labeled domains have been extensively studied in reconstructing surfaces or volumes from sets of planar cross-sections. The authors proposed a phase field framework for 3D surface reconstruction from parallel slices based on the modified Cahn-Hilliard equation [3][4][5]. All of the abovementioned reconstruction methods can process one region or separated regions but cannot process contacting regions. If we reconstruct surfaces from each single-labeled point, the resulted objects may overlap with each other (see Figure 1). Here, the surface/volumesoverlapping problem refers to the phenomenon that the volume voxel is defined for more than one value for the label. Hence, we need to develop an efficient algorithm to reconstruct multicomponent surfaces from point clouds without overlapping regions. In this paper, we propose an efficient multi-reconstruction method from point clouds. The proposed method is based on the Allen-Cahn (AC) equation [7], which has the motion of the mean curvature flow [8][9][10][11]. Under the mean curvature flow, the reconstructed surface is smooth. The Allen-Cahn equation has been applied to a wide range of problems such as crystal growth [12][13][14], image analysis [15][16][17][18], triply-periodic minimal surface [19,20], and volume reconstruction [21][22][23], and topology optimization [24]. In this paper, we use the Allen-Cahn equation in different contexts and applications. To our best knowledge, the presented approach is the first algorithm using the modified vector-valued Allen-Cahn equation for multi-reconstruction from point clouds. Our method can guarantee the production of well-reconstructed multicomponent surfaces without self-intersections. Our main contributions in this paper are as follows: • The first algorithm using a modified vector-valued Allen-Cahn equation for multireconstruction from point clouds, which can reconstruct multicomponent surfaces without self-intersections; • Based on operator splitting techniques, the proposed numerical scheme is simple to implement; • The algorithm can be straightforwardly applied to a graphics processing unit (GPU), allowing for accelerated implementation that performs many times faster than other central processing units (CPU).
The remainder of the paper is organized as follows. Section 2 briefly reviews some related works. In Section 3, we introduce the preliminary. In Section 4, the governing equations for the multi-reconstruction from a point cloud are presented. Section 5 presents the operator splitting scheme, and Section 6 states the numerical results of several computational examples. Conclusions are drawn in Section 7.

Implicit Surface Reconstruction
Implicit methods usually deal with data imperfections [25,26]. Most of the implicit implemented function are based on the signed distance function, radial basis functions, piecewise polynomial functions, or indicator functions. Hoppe et al. [27] defined the signed distance field by using the tangent plane projection. In [28], Carry et al. employed the polyharmonic radial basis function to reconstruct the manifold surfaces. This technique is well-suited to deal with non-uniformly sampled point clouds. However, it fails to reconstruct a proper surface from the point clouds with noises. Ohtake et al. [29] used the piecewise quadratic function and partition-of-unity weights to capture the local surface shape. This approach has the benefit of sharp feature reconstruction. Manson et al. [30] presented a streaming method for surface reconstruction from large point sets by using the wavelets method. The Poisson method [6] is a popular implicit method because of its resilience to data noises. However, it tends to smooth geometric details. To overcome this defect, Kazhdan and Hoppe [31] extended the Poisson method by explicitly incorporating the points as interpolation constraints by formulating a screened Poisson surface reconstruction method. Bretin et al. proposed a variational approach for the reconstruction of a volume from slices by using a minimizer of a geometric regularity criterion, Willmore energy, with inclusion-exclusion constraints associated with the cross-sections. In [32], the authors considered the problem of 3D shape reconstruction from multimodal data, given uncertain calibration parameters. To analytically and compactly represent the object, they exploited a parametric level set method, which utilized ellipsoidal radial basis functions. Kim and Lee proposed a modified Cahn-Hilliard equation for three-dimensional volume reconstruction from 2D slices [4]. In order to accurately satisfy the constraints while obtaining a smooth result, they applied a presmoothing procedure based on anisotropic diffusion to the slices.

Multi-Reconstruction
Volume fractions and implicit functions are typical approaches for defining the multilabeled domains. Hirt and Nichols [33] developed an efficient interface reconstruction method by treating complicated free boundary configurations as the material interfaces. To obtain the finer interface details, the authors utilized information about the position of the material centroid for the interface reconstruction [34][35][36]. Lemoine et al. [37] proposed a fast analytic reconstruction formula replacing the original minimization stage in the framework of the moment-of-fluid method. This algorithm can provide accurate results with a lower computational cost. Kikinzon et al. [38] extended the moment-of-fluid method to simultaneously establish the geometry and topology. Since the modified moment-of-fluid method does not require data from its neighbors, the modified moment-of-fluid method is well suitable for massively parallel computing. Yuan et al. [39] introduced a class of objectspace multiphase implicit functions that are capable of accurately modeling objects with multiple internal regions. Zhang and Qian [40] resolved the topology ambiguity method by analyzing a hybrid Octree method with both cubic and tetrahedral leaf cells. This method can apply to both homogeneous and multimaterial domains. Da et al. [41] presented the first non-manifold triangle mesh tracking method to simultaneously maintain intersectionfree meshes. Their methods support multimaterial remeshing and topological operations. In the context of multi-reconstruction from cross-sections, Liu et al. [42] proposed an extension of the projection-based approach that can deal with curve networks of arbitrary topology on cross-section planes with arbitrary orientations. Bermano et al. [43] used multiple implicit functions to extract multi-labeled material interfaces from sampled planar cross-sections of arbitrary orientation. Huang et al. [44] developed the first reconstruction algorithm for multi-labeled material interfaces that allows for explicit topology control. Li et al. [45] constructed 3D models of multi-labeled objects from sets of cross-sections by employing a multicomponent Cahn-Hilliard system. This algorithm is simple to implement and produces smooth surfaces. However, all the mentioned multi-labeled material reconstruction methods can not directly reconstruct surfaces from multi-labeled point clouds. To our best knowledge, the presented approach is the first algorithm using the modified vector-valued Allen-Cahn equation for multi-reconstruction from point clouds. Our method can guarantee the production of well-reconstructed multicomponent surfaces without self-intersections.

Preliminary: Poisson Surface Reconstruction
The Poisson reconstruction method creates smooth surfaces [6] by using an indicator ψ (defined as 1 at the voxel inside the solid M and 0 at the voxel outside the solid M).
Given the vector field v(x) ∈ R 3 defined by the data samples S, the goal of the Poisson reconstruction method is to find the scalar function ψ(x) by minimizing Here, x ∈ Ω and Ω ⊂ R 3 . ∇ is the gradient operator and · 2 is the L 2 inner product. By applying the divergence operator to form the Poisson equation, the governing equation can be obtained by where ∆ and ∇· are the Laplacian and divergence operators, respectively. The vector field v is defined by convolving the normal field with a Dirac delta function: where n(p) is the inward surface normal at p ∈ S, δ is the Dirac delta function, and ∂M is the boundary of the solid M. After solving Equation (2), we can obtain the solution ψ, which is also called the indicator function [46][47][48][49][50]. Then, we define the reconstructed surface S as the half isosurface of the scalar function ψ.

Methodology
Using the Poisson surface reconstruction [6], the indicator function of ith material can be determined by a discrete function ψ i (x) in the domain Ω. Here, ψ i (x) = 1 if the voxel x is determined as a volume voxel, otherwise ψ i (x) = 0. As mentioned in the above section, ψ i contains lots of good voxels and few self-intersection regions. We want to obtain a new smooth discrete function φ i (x), which approaches the given ψ i (x) without self-intersections. To realize this goal, we propose the following total free energy Here, the double-well potential function F(φ i ) = φ 2 i (φ i − 1) 2 /4 can enforce φ to be approximately 0 or 1. is the positive constant and relates to the phase transition width. λ is the positive constant. Minimizing the first two terms of the functional E (φ) is equivalent to finding the steady-state solution of the gradient flow equation: Here, t is the time variable. Equation (5) is the well-known Allen-Cahn equation [7], which has the motion of mean curvature [8]. This ensures that the obtained surface is smooth. Minimizing the third term of the functional E (φ) enforces that the new smooth discrete function φ i (x) is almost the same as those in the originally known volume ψ i .
The sum of the different material components should be equal to 1, i.e., ∑ N i=1 φ i = 1, which implies that any voxel is defined as one label (See Figure 2). Therefore, if we minimize the functional E and keep ∑ N i=1 φ i = 1, the obtained surface will be similar to ψ i and will be smooth without self-intersections. To minimize the proposed functional E , we consider the gradient descent method in the L 2 space: where This Lagrangian multiplier guarantees that the constraint The total energy can be decreased over time as follows: Here, we have used the periodic boundary conditions for φ i . It should be noted that the solution of Equations (6)- (8), which are derived from a constrained gradient flow of the free energy functional, is unique.

Numerical Method
We utilize an efficient operator splitting scheme for Equation (6) on the discrete domain Ω = (0, L x ) × 0, L y × (0, L z ). Let h = L x /N x = L y /N y = L z /N z be the uniform grid size. Let x p = ph, y q = qh, z r = rh, 1 ≤ p ≤ N x , 1 ≤ q ≤ N y , and 1 ≤ r ≤ N z . Let x pqr = (x p , y q , z r ) and φ k i,pqr be approximations of φ i (x pqr , k∆t), where ∆t is the time step. We split the original problem (6) into a sequence of simpler problems as Here, φ 1 i , φ 2 i , and φ 3 i denote the solutions of the subproblem (11)-(13), respectively. For a fixed x, we can see that Equation (11) is a separable ordinary differential equation, i.e., λdt + dφ/(φ − ψ) = 0. We can obtain the following solution: Next, we employ an implicit method to Equation (12): This discrete scheme can be solved by a fast Fourier transform as follows. The discrete Fourier transformφ i,lmn for l = 1, . . . , N x , m = 1, . . . , N y , and n = 1, . . . , N z is defined aŝ where ξ l = 2π(l − 1)/L x , η m = 2π(m − 1)/L y , and γ n = 2π(n − 1)/L z , respectively. The inverse discrete Fourier transform is We apply the Fourier transform to both sides of Equation (15): Then, we can get For a fixed x, Equation (13) is also a separable ordinary differential equation, i.e., With the initial condition φ 2,k+1 pqr , we can obtain the following solution after ∆t: Then, we set φ k+1 i = φ 3,k+1 i and complete the one time step processing. Finally, the proposed scheme can be written as follows: Therefore, the numerical scheme (Equations (21)- (23)) is simple to implement because Equation (6) is divided into one linear equation and two nonlinear equations. The linear equation is discretized using an implicit method, and the resulting discrete system of equations is solved by a fast Fourier transform. The two nonlinear equations are solved analytically due to the availability of a closed-form solution. Our method is fast since the computational complexities of two nonlinear equations are O(N x N y N z ) and the computational complexity of the linear equation is O(N x N y N z log (N x N y N z )). Therefore, our method is simple and fast. Our method also has merit in that it can be straightforwardly applied to GPU-accelerated implementation by using the MATLAB Parallel Computing Toolbox. The operator splitting-based hybrid numerical scheme (Equations (21)-(23)) has a good numerical stability. For Equation (21), we get 0 ≤ φ 1,k+1 (22) is a heat-type equation, its implicit numerical scheme can allow the use of a large time step. For Equation (23), we get Thus, if φ 2,k+1 i ∈ [0, 1], then φ k+1 i ∈ [0, 1]. Therefore, the numerical solution φ i is bounded, which implies that our proposed scheme allows using a larger time step. The hybrid numerical scheme (Equations (21)-(23)) has first-order accuracy in time and second-order accuracy in space because Equation (22) is employed by using an implicit method and Equations (21) and (23) are two ordinary equations. It should be noted that in the analysis, we can not prove the unconditional stability of our method. On the other hand, for the computational simulation, the results with a large time step are generally less accurate than those obtained by using a small time step, because using a larger time step would cause a large error of the numerical solutions. Therefore, a small time step will be used for highly accurate numerical solutions. Therefore, to maintain the numerical accuracy and reduce computational costs, we propose to use ∆t = 0.25 in this paper.

Experimental Results
In this section, various numerical results will be performed to demonstrate the accuracy and efficiency of our algorithm. We consider the numerical result as a steadystate solution if the relative error for every component is less than a tolerance tol, i.e., for i = 1, . . . , N. Unless otherwise specified, we use the following parameters: h = 1, ∆t = 0.25, tol = 1e − 5, λ = 0.1, and = 0.75. The reconstructed surface S is defined as the half isosurface of the phase field φ i . Note that the original input points cloud is taken from [51].

Basic Mechanism of the Algorithm
We start with an example that illustrates the basic mechanism of our algorithm. As is shown in Figure 3a, let us consider two-labeled oriented point samples of two intersecting spheres. The radii of the red sphere and blue sphere are 100 and 75, respectively. Figure 3b,c shows the reconstructed surface obtained by the Poisson reconstruction and our method, respectively. The overlapping region exists in the result of Poisson reconstruction, as shown in the left figure of Figure 4. On the other hand, our proposed method can overcome this difficulty, as shown in the right figure of Figure 4. Figure 4 illustrates the evolution of our modified multicomponent AC equation in the cross-view. From left to right, the illustration represents t = 0, 1, and 5.5, respectively. It is obvious that our method can reconstruct multicomponent surfaces without overlapping regions.

Performance on the Complex Reconstruction
We use the proposed method to reconstruct the surface from the complex multilabeled point clouds. We compared our results with the results attained by Poisson surface reconstruction in Figure 5a,c, respectively. The overlapping regions are obvious. On the other hand, our method can avoid self-intersections and provide a high-quality repairing result, as shown in Figure 5b,d. The mesh grid size is 298 × 222 × 156. The computation of our method takes about 2.5 s to perform on the GPU.

Multicomponent Surface Reconstruction
Our method can process the multicomponent surface reconstruction from the multilabeled point clouds. Figure 6a-c shows the results with two-component, three-component, and four-component reconstructions, respectively. Here, a mesh grid 298 × 156 × 212 is used. Our proposed method achieves fast convergence with few iterations. These simulations demonstrate that our proposed method can perform well for multicomponent surface reconstruction.   Figure 8 shows our multi-reconstruction results for the happy buddha surface with different densities of input point data. Figure 8a is the initial input data. Note that for the purposes of better visualization, we displayed the points more sparsely than the original density. Observing these results in Figure 8, we can see that with a low density of input point data, our proposed method can still reconstruct the surface with finer details and features. The agreement between the results obtained by different sampling densities suggests that our method can successfully reconstruct the surface with the low sampling density.

Parameter Sensitivity Analysis
In this section, we will analyze the parameter sensitivity of our proposed algorithm. The fidelity term in Equation (6)  If λ = 0, (6) becomes the classical Allen-Cahn equation. Consequently, the interface of phase field φ i will shrink with time. That is, if λ is too small, the detailed information of the original object will be lost, as shown in Figure 9a. On the other hand, if λ is too large, the fidelity term is dominant so that φ i ≈ ψ i for any time. Thus, the reconstructed surfaces will overlap with each other, as shown in Figure 9c. Therefore, we should choose an appropriate λ. In this paper, we propose to use λ = 0.1. The role of is related to the material transition width. We take the same initial condition except for different values of . Observing the results of Figure 10, we can observe that if is small, the surface of objects is sharp and the volumes have many holes. If is too large, the surface becomes over-smooth and the volumes overlap with each other. In this paper, we propose using = 0.75.  Table 1 presents the information on the grid size, iteration numbers, CPU and GPU times. CPU(pro) and GPU(pro) are the times taken to process the multi-reconstruction by using the computer's central processing unit and graphics processing unit, respectively. The CPU and GPU times in seconds of our calculations, which are performed in MATLAB, are measured on a desktop with an NVIDIA GeForce GTX 1050Ti and a quad-core Intel Core i5. We can observe that our proposed method obtains fast convergence after a few iterations because the computational complexities of the two nonlinear equations are O(N x N y N z ) and the computational complexity of the linear equation is O(N x N y N z log (N x N y N z )). Furthermore, we can employ a GPU-accelerated method that performs many times faster than CPU-only alternatives. It should be noted that the gain is related to the speed of the computation and our method may perform well in a high-performance computing cluster system. Table 1. List of grid size, iterations, CPU times (second), GPU times (second), and errors. CPU(pro) and GPU(pro) are the times taken to process the multi-reconstruction by using the computer's central processing unit and graphics processing unit, respectively.

Comparisons with Related Works and Accuracy Test
To evaluate the accuracy of our proposed method, we compare our results with the exact solution. Here, the error for every component is defined as is the result of Poisson reconstruction from the all of the point clouds. We can see that the reconstructed surfaces are close to the exact solution, as shown in Table 2. Table 2. Accurate evaluation of our method.  Based on the Cahn-Hilliard system [18,52,53], Li et al. proposed an efficient and robust algorithm to reconstruct the volumes of multi-labeled objects from sets of cross-sections without overlapping regions, artificial gaps, or mismatched interfaces [45]. Our algorithm can also handle cross-sections wherein different regions have different labels. Firstly, we obtain the point clouds from the given set of slice data by using the MATLAB contour file (See Figure 11a,b). Then, we perform our method to reconstruct the surface, as shown in Figure 11b-d. Figure 11e-h are the two-dimensional results at the slice data z = 20, 50, 80, and 110, respectively. Here, a mesh grid 128 × 128 × 124 is used. A similar test has been performed in [45]. Observing these numerical results, we can see that our proposed method is qualitatively in good agreement with the results in [45]. Our method is simpler and easy to perform compared to Li et al.'s method because their method is a fourth-order nonlinear partial differential equation and is hard to solve in an explicit manner. On the other hand, our method is a second-order partial differential equation and is easy to solve in the splitting manner. Note that this comparison is unfair because Li et al.'s method [45] can directly reconstruct the volumes of multi-labeled objects from sets of cross-sections. However, our approach should firstly translate the slice data into the point clouds and then perform the method to reconstruct surfaces, which may lead to missing data during data transforming. To our best knowledge, the presented approach is the first algorithm using the modified vector-valued Allen-Cahn equation for multi-reconstruction from point clouds.

Conclusions and Future Work
We have shown that surface reconstruction can be expressed as a modified vectorvalued Allen-Cahn equation, which obtains a smooth surface due to the motion of the mean curvature flow, and we have demonstrated that this approach can robustly reconstruct multicomponent surfaces without self-intersections. Based on operator splitting techniques, our numerical scheme is simple to implement and achieves fast convergence, which can be regarded as an advantage for real-time application. The relative computational cost is also presented. Our algorithm can also handle cross-sections wherein different regions have different labels. Compared to Li et al.'s method [45], our method is simpler and easier to perform. This comparison is considered unfair because Li et al.'s method [45] reconstructs the multi-labeled volumes from sets of cross-sections, and our method is a multi-reconstruction from point clouds. One limitation of our proposed method is that it suffers from high memory usage. As an example, Figure 7 shows the reconstruction of a stele, constructing the voxel grid at a resolution of 642 × 476 × 484. The reconstruction was computed with 1.72 GB of RAM. Trying to compute an equivalent reconstruction with methods, such as the Octree-based approach, would require constructing at a depth of nine and may require an excess of 0.178 GB of memory. In future work, we will present a parallel Octree-based mesh refinement solver for the current algorithm.