An Improved Poisson Surface Reconstruction Algorithm based on the Boundary Constraints

—The usage of the point cloud surface reconstruction to generate high-precision 3D models has been widely applied in various fields. In order to deal with the problems of insufficient accuracy, pseudo-surfaces and high time cost caused by the traditional surface reconstruction algorithms of the point cloud data, this paper proposes an improved Poisson surface reconstruction algorithm based on the boundary constraints. For large density point cloud data obtained from 3D laser scanning, the proposed method firstly uses an octree instead of the KD-tree to search the near neighborhood; then, it uses the Open Multi-Processing (OpenMP) to accelerate the normal estimation based on the moving least squares algorithm; meanwhile, the least-cost spanning tree is employed to adjust the consistency of the normal direction; and finally a screened Poisson algorithm with the Neumann boundary constraints is proposed to reconstruct the point cloud. Compared with the traditional methods, the experiments on three open datasets demonstrated that the proposed method can effectively reduce the generation of pseudo-surfaces. The reconstruction time of the proposed algorithm is about 16% shorter than that of the traditional Poisson reconstruction algorithm, and produce better reconstruction results in the term of quantitative analysis and visual comparison.


INTRODUCTION
With the popularity of 3D reconstruction technology, the surface reconstruction has been widely applied in various fields, such as the mapping [1], driverless [2], medical technology [3], robotics [4]. For the methods of reconstructing surfaces based on the 3D point cloud data, they can be generally divided into two schemes: the local methods and the global methods. The local reconstruction methods divide all the point clouds into small blocks of data, which are then reconstructed locally and finally connect them together by using some kinds of stitching functions [5,6]. This scheme retains the surface texture features well, but it is heavily sensitive to the noise. The global surface reconstruction is an approximation of all the point clouds, and the implicit surface is optimally reconstructed by solving extreme values and other methods [7]. It is under a high overall smoothness and is suitable for the interpolation and the whole repair of irregular, non-uniform scattered data. As an implicit function-based surface reconstruction method, Poisson reconstruction method reconstructs a model with watertight closure features, geometric surface properties and detail characteristics [8]. After normal estimation on the input point cloud data, the objective of this algorithm is to estimate the indicator function of the model and extract the isosurface, and then complete the surface reconstruction using the moving cube algorithm based on the indicator function and isosurface [9]. Compared to the other algorithms, it allows a hierarchy of local basis functions to be divided and the reconstruction is projected to be a Poisson space problem, combining both of the advantages of the global and the local methods.
Recently, a large number of good works on the Poisson reconstruction has been reported by the famous scholars. For example, Kazhdan et al. proposed a screened Poisson surface reconstruction algorithm to keep the sparse point structure, and solved it by using a multiple mesh algorithm [10]. Z. Xu et al. proposed an adaptive bandwidth Gaussian kernel density estimator that facilitates the removal of the noise and anomalies in the 3D reconstruction process [11]. F. Gao et al. proposed a Poisson reconstruction algorithm based on the improved isosurface extraction, which can effectively eliminate the problems of surface holes and disconnect the surface features [12]. B. Ummenhofer et al. proposed a generalized convolutional kernel for 3D reconstruction with ConvNets from point clouds [13]. Though the above reconstruction methods have made great breakthroughs in this field, there still are some drawbacks to overcome. For example, the traditional Poisson reconstruction is prone to pseudo surface and normal inconsistencies. The time complexity of reconstruction method for improved isosurface extraction is too high. The method based on an adaptive bandwidth Gaussian kernel density estimator is not effective for the reconstruction of point cloud data with small density variation. Considering the entire above problem to deal with, this paper proposes an improved Poisson reconstruction algorithm and reduces the reconstruction time of the model to improve the reconstruction efficiency.
226 | P a g e www.ijacsa.thesai.org The remainder of this paper is organized as follows. Section II presents the proposed method of reconstruction for 3D point cloud. In Section III, the experiments are implemented to evaluate the method. Conclusions are given in Section IV.

II. WHOLE METHOD
In Fig. 1, the proposed improved reconstruction algorithm consists the following steps: the set of gradients of the indicator function is determined by the normal vector of the point cloud data, so the accuracy of the normal estimation has a significant impact on the surface reconstruction results. Compared to the traditional Poisson algorithm which mostly uses principal component analysis for normal estimation, this paper proposes an improved method. Firstly, an octree is used instead of a KDtree to search the nearest neighborhood; then the normal of the point cloud is estimated by moving least squares and accelerated by OpenMP [14], and then the normal direction is adjusted consistently by a least-cost spanning tree. The traditional Poisson reconstruction algorithm is prone to generate pseudo-surfaces. In this paper, the screened Poisson algorithm is implemented by introducing constraints on the location and gradient of the points to constrain the surface reconstruction process; and adding Neumann boundary constraints to make the solution of the indicator function more accurate and generate more accurate surface models. The whole flow of the specific algorithm is shown in Fig. 1.

A. Octree and the Moving Least Squares
The octree is used to represent a three-dimensional space, and is a spatial extension of the quadtree [15]. The geometric entities in the three-dimensional space are first dissected into cubes, each with the same time and space complexity, and then the geometric objects in the three-dimensional space of size (2n*2n*2n) are dissected by a circular recursive partitioning method to construct a directional graph with a root node [16]. In the octree structure if the cubes being divided have the same properties, the cube forms a leaf node [17]; otherwise the dissection of the cube into eight sub-cubes continues [18]. For (2n*2n*2n) size space objects are dissected at most n times, and the structure is shown in Fig. 2 [19]. The octree algorithm is more efficient than the KD-tree in searching point cloud data with high data volumes and is more automated, allowing optimization of the time for 3D reconstruction of point clouds [20,21].
The Moving Least Squares algorithm is simple and easy to implement as a method for interpolating discrete data, with high fitting accuracy [22]. When a large amount of discrete data is distributed in a heterogeneous manner, the use of traditional least squares algorithms often requires the fitting of segments to the data, in addition to avoiding the problem of discontinuous and unsmooth fitting curves on adjacent segments [23]. Moving least squares introduces the concept of tight branches and the fitting function is constructed differently.
1) The fitting function: The fitting function [24] for the moving least squares in the local area is shown in Eq. (1). In this paper, the quadratic basis is used to improve the accuracy of surface fitting and normal estimation [25]. The squared weight of the difference between the local approximation of ( ) and the nodal value is minimized, the weighted residual equation is shown in Eq. (2).
where is the number of nodes in the area of influence; ( ) is the weight function of node . In order to determine the coefficient ( ) , it must be made to take a minimal value. The derived equation is shown in Eq. (3). where The coefficients of the fitting equation are shown in the following equation.
where ( ) is the shape function; is the order of the basis function.
2) The selection of the weight function: The weight function ( ) in Eq. (2) is compactly supported in the moving least squares algorithm, that is, the weight function is only affected by the subdomain near , which is called the influence region of point . Beyond this region, its weight is small and the influence can be ignored. As an indispensable part of the moving least square method, the selection of weight function is very important to the fitting accuracy [26,27]. Generally speaking, the circle is chosen as the support domain of the weight function, and its radius is denoted as , as shown in Fig. 3. The common weight functions include the spline weight function and the Gaussian weight function [28,29]. When Gaussian weight function is used for the moving least squares surface fitting, the width of a single kernel is difficult to meet the feature requirements of the whole model. It needs to be selected according to the point cloud density and curvature distribution in the local neighborhood, and the kernel width is not easy to be determined. As reported in literature [30], the cubic spline weight function has continuity. Therefore, the cubic spline weight function is selected in this paper. Denoting , ̅ , the cubic spline weight function is shown as follows:

B. The Normal Estimation and Direction Consistency
Define a local quadratic surface as ( ) ( ( )). The first order partial derivative of the surface at a point is then calculated as shown in the following equation.
where ( ) . So the normal direction of the surface at this point is calculated by the following formula. Where, , , , , and f are the coefficients of the surface equation.
The normal direction calculated based on the moving least square method is ambiguous, that is, only the line where the normal line is located is obtained, but the direction of the normal line is not determined. Therefore, it is necessary to carry out direction consistency processing for the normal in the model to ensure its accuracy [31]. In this paper, the minimum spanning tree (MST) is used to adjust the normal direction uniformly [32,33]. Firstly, the algorithm defines a cost function for the point cloud in the model, as shown in the following equation.
Where is the unit vector pointing from point to point . and are adjacent points, and and are normal to points and , respectively.

C. Screening Factors
Because Poisson equation is a partial differential equation, some errors will cause the indicator function to shift, resulting in a large number of pseudo-closed surfaces at the edges of non-watertight surface connections [34,35]. In this paper, the process of surface reconstruction is constrained by introducing the position and gradient constraints, and then the screened Poisson reconstruction algorithm is implemented. www.ijacsa.thesai.org In the process of solving the Poisson equation, the indicator function is solved for the directed point set ⃗ , so that the gradient of the function and the directed point set form the best approximation, that is, the scale function is minimized to solve the problem ‖ ⃗ ⃗ ‖ , ( ) ∫‖ ( ) ⃗ ( )‖ . In this paper, a point set with weight is given, and some gradient constraints and value constraints of discrete points are added to minimize the scale function to recalculate the sample point function. The constraint equation is shown in Eq. (10).
where, is the sample point of the input point set, is the gradient of the indicator function, ( ) is the surface region to be reconstructed, is the screening factor, which is used to measure the fitting gradient and the fitting value, and ( ) is the sample point weight. For the convenience of calculation, the weight ( ) of each sample point is set as 1 in this paper. Eq. (11) can be obtained by simplified to be: where, represents the standard inner product on the space of functions in the unit cube (scalar and vector values). ( ) ( ) is the representation of the unit cube on the function space. This value is obtained from the weighted summation of the functions of the sampling points, as shown in Equation (12).

D. The Neumann Boundary Constraint
In order to satisfy the solution of the shielded Poisson equation and the boundary conditions in a given region, the algorithm bias error is reduced by introducing boundary constraints for solving the partial differential equation. There are three main types of boundary constraints: the Dirichlet boundary constraint [36], the Neumann boundary constraint and the Robin boundary constraint [37,38]. At the endpoints it is generally written in the form Ay + By' = C. a) If B = 0 and A ≠ 0, it is called the Dirichlet boundary condition, also called the first type of boundary condition, and gives the corresponding value of the unknown function on the boundary. b) If B ≠ 0, A = 0, then it is called Neumann boundary condition, also called the second type of boundary condition, and gives the directional derivative of the unknown function normal to the boundary outside. c) If B ≠ 0, A ≠ 0, then it is called Robin boundary condition, also called the third type of boundary condition, and gives a linear combination of the value of the function of the unknown function on the boundary and the directional derivative of the outer normal.
The Neumann boundary constraint is added to the Poisson reconstruction process to force the normal derivative of the implicit function to be zero at the boundary, which makes the solution of the indicator function more accurate. Both the Dirichlet and Neumann boundary constraints require the implicit function to take a negative value outside the model. Since the gradient of the indicator function is equal to the vector field, Neumann boundary conditions are more conducive to the solution of the indicator function than Dirichlet. Neumann boundary conditions constrain the watertight surface and require the derivative of the implicit function to be zero on the boundary of the integral domain. This property is compatible with the gradient of the indicator function, since the vector field V is numerically set to zero in the region far from the point cloud sample. When the surface exceeds the boundary of the integral domain, the Neumann boundary constraint creates a deviation across the boundary of the domain.

III. EXPERIMENTAL RESULTS AND DISCUSSION
To verify the effectiveness, accuracy and feasibility of the proposed algorithm, it is compared with the traditional Poisson reconstruction algorithm and greedy projection triangulation algorithm [39,40]. The 3D scanning datasets from Stanford University and GeometryHub are used in this paper. All experiments are performed on a computer with a quad-core AMD Ryzen 3 3100 CPU and 16GB RAM configured at 3.59 GHz. And the parallel acceleration is implemented by the OpenMP.

A. Results
In Fig. 4, the original point cloud data of Model 1, Model 2 and Model 3 are shown from left to right. Fig. 5 shows the normal estimation of the MLS algorithm for the three data respectively, and the calculated normal is ambiguous. Fig. 6 shows the initial normal is estimated by the moving least squares algorithm for three data, and then the minimum spanning tree is used to redirect the normal. Based on the point cloud data in Fig. 6, the proposed algorithm is compared with the traditional Poisson reconstruction algorithm, the greedy projection triangulation algorithm.  229 | P a g e www.ijacsa.thesai.org  Fig. 7, Fig. 7(a), 7(d) and 7(g) show the model generated by the traditional Poisson reconstruction algorithm, Fig. 7(b), 7(e) and 7(h) show the model generated by the greedy projection triangulation algorithm, and Fig. 7(c), 7(f) and 7(i) show the model generated by the proposed algorithm in this paper.
In order to further optimize the reconstruction effect, the different parameters of the proposed algorithm are analyzed and compared with the point cloud data of the above models. The effect diagram of model 1 is shown in Fig. 8. Then the above point cloud data were reconstructed with three types of boundary conditions respectively. The experimental results are shown in Fig. 9 and Table II. Fig. 9(a), 9(d) and 9(g) is the surface reconstruction under Neumann boundary constraints, while Fig. 9(b), 9(e) and 9(h) is the surface reconstruction under Dirichlet boundary constraints. Fig. 9(c), 9(f) and 9(i) shows the surface reconstruction under Robin boundary constraints.

B. Discussions
It can be seen from Fig. 5, the normal direction is confused, someone pointing to the inner side, and someone pointing to the outer side. If there is no redirection, further processing will lead to many reconstruction errors. Fig. 6 shows the use of normal redirection optimizes the consistency of the normal direction and provides a more accurate normal input for later reconstruction. As can be seen from the comparison of several groups of models in Fig. 7, the surface features of the traditional Poisson algorithm have poor sealing property, and it is easy to generate pseudo surfaces by misconnecting the regions that do not belong to morphological features at the edges, among which Fig. 7(g) is the most obvious. The model reconstructed by greedy projection triangulation algorithm has many small holes, and the model surface is rough. Compared with the former, the proposed algorithm can effectively solve the whole problem by introducing shielding factor and boundary constraint. And the overall model is more perfect in reduction and surface smoothness.
Based on the data in Table I, the reconstruction time of the traditional Poisson reconstruction algorithm, greedy projection triangulation algorithm and the proposed algorithm is compared. The reconstruction time of the proposed algorithm is about 16% shorter than that of the traditional Poisson reconstruction algorithm, but the reconstruction time is longer than that of the greedy projection triangulation algorithm. This is because the proposed algorithm uses the moving cube algorithm to extract the isosurface, which consumes a long time and needs further improvement. When calculating the point cloud model with larger data volume, the algorithm in this paper has higher reconstruction efficiency. In Fig. 8, the higher the value of the screening factor, the more detailed the reconstructed model will be, but the longer it will take. At the screening factor of 4, the detailed features are well represented and take a moderate amount of time. As can be seen from Fig. 9, the accuracy of the reconstruction model under Dirichlet boundary constraint is low, and many detailed features are lost. Although Robin boundary constraint is a combination of Dirichlet boundary constraint and Neumann boundary constraint, it takes the longest time and has a poor reconstruction effect on the edge, as shown in the Fig. 9(c) for the pseudo-surface marked by the red circle. The reconstruction time under the constraint of Neumann boundary is the shortest. Compared with the traditional Poisson surface reconstruction, the pseudo-surface is significantly reduced, and the effect is more ideal.   In this paper, we propose an improved Poisson surface reconstruction algorithm based on the Boundary Constraints. Firstly, octree is used to replace KD-tree for nearest neighbor search. Secondly, the normal vector is estimated by moving least square method, and the redirection based on minimum cost spanning tree is used to reduce the error. Finally, on the basis of traditional Poisson reconstruction, screening factor and Neumann boundary constraint are introduced to further improve the reconstruction effect.
The experimental results on different data show that the proposed algorithm achieves more accurate reconstruction results, which can effectively reduce the generation of pseudosurfaces and also reduce the running time to a certain extent. The further work is to improve the extraction of isosurface on the basis of the proposed algorithm, and try to apply the modified algorithm to other fields to obtain high-quality reconstructed models.

ACKNOWLEDGMENT
This study was supported by A project ZR2021MF017 supported by Shandong Provincial Natural Science Foundation; A project ZR2020MF147 supported by Shandong Provincial Natural Science Foundation; A project 2020SNPT0055 supported by SDUT&Zibo City Integration Development Project; The National Natural Science Foundation of China (61502282).