Fast Implicit Surface Reconstruction for the Radial Basis Functions Interpolant

: In this paper we improve an e ﬃ cient implicit surface reconstruction method based on the surface following method for the radial basis functions interpolant. The method balances the reconstruction e ﬃ ciency and the evaluation e ﬃ ciency in the process of surface following. The growing strategy of the surface following method combines both the evaluation and reconstruction processes. Based on the analysis of the black-box fast multipole method (FMM) operations, we improve the FMM procedures for single point evaluation. The goal is to ensure that one point evaluation of the method obtains an optimum e ﬃ ciency, so that it can be e ﬃ ciently applied to the voxel growing method. Combined with the single point FMM, we improve the voxel growing method without manually specifying the seed points, and the leaf growing method is developed to avoid a mass of redundant computation. It ensures a smaller number of evaluation points and a higher evaluation e ﬃ ciency in surface following. The numerical results of several data sets showed the reliability and performance of the e ﬃ cient implicit surface reconstruction method. Compared with the existing methods, the improved method performs a better time and space e ﬃ ciency.


Introduction
The reconstruction of implicit surfaces is commonly used in computer graphics, scientific data visualization, geometry processing, reverse engineering, etc. The reconstruction efficiency has always been a difficult problem that hinders the wide application of implicit surface (e.g., real-time shape editing). This work tries to extract the isosurface from the zero level set of an implicit function efficiently.
As an important interpolation method, the radial basis functions (RBFs) interpolant [1][2][3] is widely used in different applications. In numerical analysis, the RBF interpolant plays a key role in solving partial differential equations (PDEs) [4]. In nonlinear time series analysis, it is used for the approximation of the flow in the phase space [5,6]. In spatial interpolation, it is used to interpolate different geological data [7] (e.g., grade, fault and formation).
In implicit modeling, given a set of unorganized points with normals (Hermite data), the implicit surface reconstruction methods interpolate the unknown surface using an implicit function. Based on the idea of signed distance field interpolation [8], the radial basis functions interpolant can be used to implicitly represent the implicit surface satisfying the Hermite data. The unorganized points are used to construct on-surface constraints, and the offset points generated by projecting along the normals are used to construct off-surface constraints.
The implicit surface reconstruction involves two main processes, namely evaluation and reconstruction. As a classical method of surface reconstruction, the marching cubes (MC) method [9,10] extracts the isosurface by triangulating the cubes in the discretized space. The method evaluates the implicit function values of all points sampled on the regular voxel points. It takes too much time to evaluate the invalid voxels that do not intersect the isosurface. The improved surface following method [11,12] minimizes the number of evaluation points by tracking voxels intersecting the isosurface.
For the RBF interpolant with large numbers of constraints, both the evaluation and reconstruction processes should be combined to accelerate the reconstruction of the implicit surface. The fast multipole method (FMM) [13,14] is performed to efficiently calculate the sums of pairwise interactions between source points. However, different from the marching cubes method, the evaluation points are determined dynamically by tracking seed voxels for the surface following method. Therefore, on the one hand, the original FMM should be improved to efficiently evaluate the points one by one. On the other hand, to avoid the repeat computation of the coefficients of the FMM, the evaluation points in the same leaf cell of the FMM tree should be evaluated first in the reconstruction process.
According to the above analysis, the minimum evaluation points and the optimization of single point evaluation should be combined to improve the efficiency of implicit surface reconstruction. In this paper, we improve an efficient implicit surface reconstruction method for the RBF interpolant. The method balances the reconstruction efficiency and the evaluation efficiency in the process of surface following. The goal is to ensure that one point evaluation of the method obtains a similar efficiency to the original FMM, so that it can be efficiently applied to the voxel growing method. To efficiently evaluate the field points one by one in reconstruction, we improve the FMM procedures for single point evaluation. The numerical results of several data sets showed the reliability and performance of the efficient implicit surface reconstruction method.

Related Works
We briefly review the two bodies of work that are close to ours, namely radial basis functions and implicit surface reconstruction.
There are many implicit functions (e.g., Poisson function [15,16], the Discrete Smooth Interpolation function [17,18]) used to represent implicit surfaces. For the radial basis functions interpolant, the evaluation process is to compute the sums of the linear combinations of kernel functions. The common method for fast evaluation is the fast multipole method (FMM). The basic idea is to expand the radial basis function (RBF) kernels via low-rank approximations, and neglect the far field expansion within a given precision. At present, a number of evaluation methods have been proposed to expand different types of kernels in an analytic way, such as thin-plate splines [19], Gaussian kernel [20] and polyharmonic splines [21]. Recently, some general FMM methods are developed to evaluate kernels without the implementation of multipole expansions, including the kernel independent FMM [22] and the black-box FMM [23]. Besides the global evaluation methods, the local evaluation methods using compactly supported radial basis functions (CSRBF) are also a feasible way of acceleration. To balance the efficiency and effect, the local evaluation methods require the selection of adaptive supported radii for different interpolation centers. Based upon the comparison of a quantitative measure of approximation error, Zhang et al. [24] proposed an adaptive radial basis functions interpolation method using an error indicator.
The reconstruction process usually has nothing to do with the evaluation process. Some of the implicit surface reconstruction methods only need to know the implicit function values at any given position, such as the well-known marching cubes and marching tetrahedra algorithms. At present, many extensions of marching cubes have been proposed, such as dual marching cubes [25], adaptive marching cubes [26], extended marching cubes [27] and multiple material marching cubes [28]. To track the zero level sets of the implicit surface directly, some surface following methods are developed based on a greedy growing idea. Lee et al. [11] presented the growing-cube algorithm to track the neighboring cells with zero values based on the marching cubes algorithm.
Note that the combination of some additional conditions (e.g., the Hermite data of the implicit function) is useful to improve the mesh quality and minimize the number of evaluations.
Treece et al. [29] presented the regularized marching tetrahedra algorithm to reconstruct meshes with fewer triangles and better aspect ratios. On purely the reconstruction process, the surface following method with mesh optimization strategy is an optimum scheme for implicit surface reconstruction. In this work, we focus on improving the reconstruction efficiency, combining both the evaluation and reconstruction processes.

Implicit Function
In implicit modeling, the implicit surface S is defined as the zero level set x f (x, y, z) = 0 of the implicit function f (x). To reconstruct the implicit surface S, the implicit function f (x) should be first obtained by interpolating the domain constraints where f i are function values of the geometry domain and N is the number of constraints. These interpolation constraints are usually constructed from a set of sampling points. For the radial basis functions interpolant, the implicit function f (x) has the form [3] f where x = (x, y, z), K x, y j is a kind of globally supported radial basis functions and p(x) are low-order polynomials. The unknown coefficients ω i can be determined by solving the linear system combined by the domain constraints. The first variable x in K(x, y) is viewed as an evaluation point (or target point), and the second variable y is viewed as a source point. Taking the polynomial part p(x) = c 1 + c 2 x + c 3 y + c 4 z as an example, the smoothest interpolant satisfies the orthogonality conditions These orthogonality conditions, along with the domain conditions, lead to a linear system. The matrix form of the linear system can be written as where K i,j = K x i , y j , and the unknown coefficients ω i and c i can be determined by solving the equation.
There are several ways (e.g., Fast RBF method [30,31]) to efficiently solve the interpolation equation for the RBF interpolant. In this paper, we focus more on the efficient reconstruction of the implicit surface. As mentioned earlier, the evaluation process and the reconstruction process are combined to improve the reconstruction efficiency.

Black-Box FMM
To efficiently evaluate the function values of the RBF implicit function f (x), the fast multipole method is used to compute the sums of RBFs based on the low-rank approximations. Taking the 1-D black-box FMM [23] as an example, the background of the kernel independent FMM is reviewed briefly.
Given a set of target points {x i } and source points y j , the kernel K(x, y) is expanded in a low-rank approximation using Chebyshev polynomials where p is the expansion order of series, x n and y m are Chebyshev nodes, S p (x n , x) is the interpolation function at node x n , and S p y m , y is the interpolation function at node y m . To make the expansions with a rapid convergence, the evaluation points are well separated from the source points (far field). If the evaluation points and source points (near field) are not well separated, the kernels will be evaluated directly. Then the sum of kernels can be efficiently evaluated by changing the order of computation where the number of source points N satisfies N = N f ar + N near . To divide the near field and far field adaptively, the multilevel FMM approach [32] is developed for hierarchical evaluation using the hierarchical data structure, as shown in Figure 1. The source points are divided into rectangle boxes in different levels using the Quadtree data structure. The well-known FMM operations of the near field and far field are divided into several steps, including P2M, M2M, M2L, L2L, L2P and P2P. The operations correspond to a tree structure.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 4 of 14 1-D black-box FMM [23] as an example, the background of the kernel independent FMM is reviewed briefly. Given a set of target points { } and source points { }, the kernel ( , ) is expanded in a low-rank approximation using Chebyshev polynomials where is the expansion order of series, { ̅ } and { ̅ } are Chebyshev nodes, ( ̅ , ) is the interpolation function at node ̅ , and ( ̅ , ) is the interpolation function at node ̅ .
To make the expansions with a rapid convergence, the evaluation points are well separated from the source points (far field). If the evaluation points and source points (near field) are not well separated, the kernels will be evaluated directly. Then the sum of kernels can be efficiently evaluated by changing the order of computation where the number of source points satisfies = + . To divide the near field and far field adaptively, the multilevel FMM approach [32] is developed for hierarchical evaluation using the hierarchical data structure, as shown in Figure 1. The source points are divided into rectangle boxes in different levels using the Quadtree data structure. The well-known FMM operations of the near field and far field are divided into several steps, including P2M, M2M, M2L, L2L, L2P and P2P. The operations correspond to a tree structure.

Single Point FMM
The original FMM evaluates the sums of matrix-vector products. For the surface following method, the evaluation points are determined via the tracking process step by step. To efficiently evaluate the field points one by one, the tree-based procedures of the multilevel approach should be improved for single point evaluation.
For a leaf cell T in the FMM tree (e.g., quad-tree in 2D, octree in 3D), the evaluation of any target point in T is decomposed into three steps based on the black-box FMM operations [23,33].

Single Point FMM
The original FMM evaluates the sums of matrix-vector products. For the surface following method, the evaluation points are determined via the tracking process step by step. To efficiently evaluate the field points one by one, the tree-based procedures of the multilevel approach should be improved for single point evaluation.
For a leaf cell T in the FMM tree (e.g., quad-tree in 2D, octree in 3D), the evaluation of any target point x i in T is decomposed into three steps based on the black-box FMM operations [23,33].
1. Set up (P2M and M2M). The P2M and M2M operations have no connection with the evaluation points and can be computed first. For any leaf cell T , the multipole moments of T are interpolated from the source points inside T at the Chebyshev nodes y T For any non-leaf cell T, the multipole moments of T are computed by gathering the child cells' multipole expansions interpolated at the Chebyshev nodes y T where T are the child cells of T and m is associated with T .
where J are the cells in the interaction list of T. The M2L operation of the computed cell is marked as true to avoid repeated calculation. For the cells T containing x i in different octree levels, the local expansions of T are computed by adding the far field interactions and translating the parent cell's local expansions interpolated at the Chebyshev nodes x T n L T n = L T,J n + n L T n S p x T n , x T n , n = 1, 2, . . . , p where T is the parent cell of T and n is associated with T . The L2L operation of the computed cell is marked as true to avoid repeated calculation. 3. Evaluation (L2P and P2P). Once the setup and translation steps have been completed, the evaluation for a single point can be computed as the sum of the far field and near field by For any target point x i in leaf cell T, the interactions of far field are computed by For any target point x i in leaf cell T, the interactions of near field are directly computed by where J are the cells in the neighbor list of T. The setup step computes the multipole moments of the cells in each level via the source points. It is only executed once before the evaluation process without considering the evaluation points. The other two steps are computed based on the evaluation points. The translation step computes the conversion of multipole moments and local expansions to the target leaf cell. Therefore, it is only executed once for the evaluation points in the same leaf cell.
For the evaluation points in different leaf cells, the computed M2L and L2L operations are also reusable without repeated calculation. The evaluation step is always executed for each evaluation point.

Voxel Growing
To polygonize the surface of an implicit function, the geometry space should be discretized to construct the voxels for evaluation, as shown in Figure 2. The classical marching cubes method evaluates all the voxels to extract the zero level set (Figure 3a), while the voxel growing method only evaluates the voxels intersecting the isosurface in the geometry space ( Figure 3b). However, it is time consuming to determine the first initial seed point by traversing all the cubes. Moreover, the original voxel growing method cannot ensure to recover the complete model with multiple isosurfaces for missing some possible seed points, as shown in Figure 3b.
To polygonize the surface of an implicit function, the geometry space should be discretized to construct the voxels for evaluation, as shown in Figure 2. The classical marching cubes method evaluates all the voxels to extract the zero level set (Figure 3a), while the voxel growing method only evaluates the voxels intersecting the isosurface in the geometry space ( Figure 3b). However, it is time consuming to determine the first initial seed point by traversing all the cubes. Moreover, the original voxel growing method cannot ensure to recover the complete model with multiple isosurfaces for missing some possible seed points, as shown in Figure 3b.  In this paper, combined with the single point FMM, we improve the voxel growing method without specifying the seed points manually. For the RBF interpolant, all of the domain constraints with zero values are selected as seed points to avoid missing evaluations, as shown in Figure 2. The To polygonize the surface of an implicit function, the geometry space should be discretized to construct the voxels for evaluation, as shown in Figure 2. The classical marching cubes method evaluates all the voxels to extract the zero level set (Figure 3a), while the voxel growing method only evaluates the voxels intersecting the isosurface in the geometry space ( Figure 3b). However, it is time consuming to determine the first initial seed point by traversing all the cubes. Moreover, the original voxel growing method cannot ensure to recover the complete model with multiple isosurfaces for missing some possible seed points, as shown in Figure 3b. In this paper, combined with the single point FMM, we improve the voxel growing method without specifying the seed points manually. For the RBF interpolant, all of the domain constraints with zero values are selected as seed points to avoid missing evaluations, as shown in Figure 2. The voxel containing at least one seed point is viewed as an initial seed voxel. Then the method tracks In this paper, combined with the single point FMM, we improve the voxel growing method without specifying the seed points manually. For the RBF interpolant, all of the domain constraints with zero values are selected as seed points to avoid missing evaluations, as shown in Figure 2. The voxel containing at least one seed point is viewed as an initial seed voxel. Then the method tracks the voxels intersecting the isosurface using all the initial seed voxels. We construct a voxel growing queue to store the adjacent voxels intersecting the isosurface for tracking. The adjacent voxels intersected the isosurface are appended to the voxel growing queue greedily until there are no new voxels intersected the isosurface, as shown in Figure 3d. Note that the evaluated voxels should be marked to avoid recomputation. Different from the marching cubes method, the voxel growing method evaluates the voxels greedily. Though the voxel growing method reduces the number of evaluations, the evaluation points should be determined via the reconstruction process step by step, which reduces the efficiency of single point evaluation. To improve the evaluation efficiency, the single point FMM is used to evaluate the points one by one. Note that the evaluated points should be marked to avoid recomputation.
The improved voxel growing achieves the same reconstruction effect more efficiently. In the growing process, a voxel growing queue is used to track an isosurface. Reconstructing an isosurface requires only one initial seed voxel. As long as each isosurface has at least one constraint point with zero values, the original model can be completely reconstructed. Generally, the sampling points can meet the requirement in practical applications. It is worth noting that the initial seed voxels in the same isosurface are repeated and should be removed to avoid recomputation. Considering the fact that a voxel growing queue corresponds to an isosurface, an available approach is to remove the repeated initial seed voxels in the same growing queue gradually.

Leaf Growing
In the voxel growing queue, the numerical results show that assigning different priorities to the pending voxels leads to a different evaluation efficiency. It is because the order of voxel growing affects the single point evaluation efficiency. Compared with the original FMM, the third step of the single point FMM computes the L2P coefficients many times for the evaluation points in the same leaf cell. For the evaluation points in the same leaf cell, the L2P operation interpolates at the same Chebyshev nodes, and the same coefficients can be precomputed by expanding the Chebyshev polynomials.
To avoid a mass of redundant computation, we present the leaf growing method based on the improved voxel growing method. The main difference is that the growing process considers both the evaluation cell of FMM and the reconstruction cell of MC. The evaluation cell of FMM refers to leaf cell, the finest level of octree. The reconstruction cell of MC refers to the voxel or cube. An available approach is to track the implicit surface leaf by leaf instead of voxel by voxel, and then to evaluate all the voxels in the same leaf.
To reduce the number of evaluation points, the improved growing strategy is to track the implicit surface leaf by leaf in global, and voxel by voxel in local. In addition, the voxels in the growing queue are specified a priority via the position in the leaf. There are three positional relationships between the leaf and the voxel, namely interior, intersection and exterior, as shown in Figure 4. In the voxel growing process, a leaf T with the precomputed L2P coefficients is viewed as a current leaf. The precomputed L2P coefficients of T will be stored until there are no voxels inside or intersecting T in the voxel growing queue. To save the memory storage, for the voxels in the growing queue, the voxels inside the current leaf should be evaluated first. Then the voxels intersecting the current leaf will be evaluated and a new current leaf will be obtained. Given a radial basis functions interpolant ( ), and a desired reconstruction accuracy , based on the above analysis, a simplified procedure of the fast implicit surface reconstruction method for the RBF interpolant is given below.
Step 1: Space division. Compute the minimum bounding box via the constraint points and divide the space into regular voxels (cubes) via the desired reconstruction accuracy . Scale the minimum bounding box in a certain ratio.
Step 2: Construct the initial seed voxels. Select the points of the domain constraints with zero values ( ( ) = 0) as the seed points and construct an initial seed voxel set via the positions of the seed points.
Step 3: Execute the setup step. Compute the multipole moments based on the source points { }. The source points are the interpolation centers for the RBF interpolant ( ).
Step 4: Traverse the initial seed voxel set . Construct a voxel growing queue and add an initial seed voxel to .
Step 5: Traversing the voxel growing queue , take out a voxel with the highest priority. Evaluate the points of the voxel by executing the translation and evaluation steps. Obtain the current leaf cell of the FMM according to the L2P operation. Store the evaluated points and the evaluated voxels to avoid the repeated evaluation.
Step 6: According to the lookup table of the marching cubes method, triangulate the voxel via the eight evaluated values. Optionally, some additional strategies can be implemented to resolve the ambiguities.
Step 7: Determine the new growing voxels for tracking. Among the six adjacent voxels of , append the un-evaluated voxels that are intersected the isosurface to the voxel growing queue .
Step 8: Specify a priority for the new growing voxels via the corresponding position in the current leaf cell of the FMM.
Step 9: Compare the new growing voxels with the initial seed voxel set , and remove the duplicate initial seed voxels in .
Step 10: Return to Step 5 until the traversing of the voxel growing queue is finished.
Step 11: Extract an isosurface based on the triangulated voxels.
Step 12: Return to Step 4 until the traversing of the initial seed voxel set is finished.
Step 13: Output the reconstructed isosurfaces. Given a radial basis functions interpolant f (x), and a desired reconstruction accuracy , based on the above analysis, a simplified procedure of the fast implicit surface reconstruction method for the RBF interpolant is given below.
Step 1: Space division. Compute the minimum bounding box via the constraint points and divide the space into regular voxels (cubes) via the desired reconstruction accuracy . Scale the minimum bounding box in a certain ratio.
Step 2: Construct the initial seed voxels. Select the points of the domain constraints with zero values ( f (x i ) = 0) as the seed points and construct an initial seed voxel set S via the positions of the seed points.
Step 3: Execute the setup step. Compute the multipole moments based on the source points y j . The source points are the interpolation centers for the RBF interpolant f (x).
Step 4: Traverse the initial seed voxel set S. Construct a voxel growing queue G and add an initial seed voxel to G.
Step 5: Traversing the voxel growing queue G, take out a voxel V with the highest priority. Evaluate the points of the voxel V by executing the translation and evaluation steps. Obtain the current leaf cell of the FMM according to the L2P operation. Store the evaluated points and the evaluated voxels to avoid the repeated evaluation.
Step 6: According to the lookup table of the marching cubes method, triangulate the voxel V via the eight evaluated values. Optionally, some additional strategies can be implemented to resolve the ambiguities.
Step 7: Determine the new growing voxels for tracking. Among the six adjacent voxels of V, append the un-evaluated voxels that are intersected the isosurface to the voxel growing queue G.
Step 8: Specify a priority for the new growing voxels via the corresponding position in the current leaf cell of the FMM.
Step 9: Compare the new growing voxels with the initial seed voxel set S, and remove the duplicate initial seed voxels in S.
Step 10: Return to Step 5 until the traversing of the voxel growing queue is finished.
Step 11: Extract an isosurface based on the triangulated voxels.
Step 12: Return to Step 4 until the traversing of the initial seed voxel set is finished.

Numerical Results
We implemented the improved implicit surface reconstruction method for the RBF interpolant and tested the performance on a Windows 64-bit PC with 3.00 GHz Intel(R) Core(TM) i5-7400 and 8GB RAM. The ScalFMM library was used to implement the single point FMM for fast evaluation.

Performance
To validate the performance of the single point FMM, a series of numerical experiments were conducted using the randomly generated source points and evaluation points in an ellipsoid. The same parameters were used for all the experimental results: the expansion order is p = 10 and the kernel is biharmonic spline. The number of evaluation points varies from 10 2 to 10 7 .

Numerical Results
We implemented the improved implicit surface reconstruction method for the RBF interpolant and tested the performance on a Windows 64-bit PC with 3.00 GHz Intel(R) Core(TM) i5-7400 and 8GB RAM. The ScalFMM library was used to implement the single point FMM for fast evaluation.

Performance
To validate the performance of the single point FMM, a series of numerical experiments were conducted using the randomly generated source points and evaluation points in an ellipsoid. The same parameters were used for all the experimental results: the expansion order is = 10 and the kernel is biharmonic spline. The number of evaluation points varies from 10 2 to 10 7 .    Appl. Sci. 2019, 9, x FOR PEER REVIEW 9 of 14

Numerical Results
We implemented the improved implicit surface reconstruction method for the RBF interpolant and tested the performance on a Windows 64-bit PC with 3.00 GHz Intel(R) Core(TM) i5-7400 and 8GB RAM. The ScalFMM library was used to implement the single point FMM for fast evaluation.

Performance
To validate the performance of the single point FMM, a series of numerical experiments were conducted using the randomly generated source points and evaluation points in an ellipsoid. The same parameters were used for all the experimental results: the expansion order is = 10 and the kernel is biharmonic spline. The number of evaluation points varies from 10 2 to 10 7 . Figures 5-7 demonstrate comparison results of the evaluation time using the original FMM (blue) and the single point FMM (red) with different steps. Note that the original FMM evaluates all the points as a whole in parallel. For the single point FMM, only the setup step is parallel. The comparison results show that the original FMM evaluating all the points as a whole has the best efficiency. If the points should be evaluated one by one, it is necessary to improve the FMM operations to achieve an optimum efficiency. The three steps of the single point FMM optimize the evaluation process by reducing the redundant computation as much as possible. The setup step improves the efficiency by precomputing the P2M and M2M operations. The translation step improves the efficiency by reducing the M2L and L2L operations, and the evaluation step improves the efficiency by optimizing the L2P operation.

Case Studies
We sampled sets of unorganized points with normals from a number of real objects to test the implicit surface reconstruction method, as shown in Figure 8. The unorganized points are used to construct the on-surface constraints ( ( ) = 0) and the normals are used to construct off-surface constraints ( ( ) > 0 or ( ) < 0) by projecting the points along the normal direction. The implicit function (the RBF interpolant) is obtained by interpolating the constraints and solving the corresponding linear system. Then the implicit surface can be reconstructed using the improved surface reconstruction method. In the process of surface reconstruction, all of the on-surface constraints are selected as the seed points. Figure 9 demonstrates comparison result of the voxel growing method with one and multiple seed points based on the on-surface constraints. The on-surface constraints are constructed by the sampling points. As mentioned, for the surface following method, the main problem is the incomplete reconstruction of sub-objects. Compared with the original voxel growing method, the improved method tends to recover the complete model as long as the sampling is sufficient.

Case Studies
We sampled sets of unorganized points with normals from a number of real objects to test the implicit surface reconstruction method, as shown in Figure 8. The unorganized points are used to construct the on-surface constraints ( f (x i ) = 0) and the normals are used to construct off-surface constraints ( f (x i ) > 0 or f (x i ) < 0) by projecting the points along the normal direction. The implicit function (the RBF interpolant) is obtained by interpolating the constraints and solving the corresponding linear system. Then the implicit surface can be reconstructed using the improved surface reconstruction method.

Case Studies
We sampled sets of unorganized points with normals from a number of real objects to test the implicit surface reconstruction method, as shown in Figure 8. The unorganized points are used to construct the on-surface constraints ( ( ) = 0) and the normals are used to construct off-surface constraints ( ( ) > 0 or ( ) < 0) by projecting the points along the normal direction. The implicit function (the RBF interpolant) is obtained by interpolating the constraints and solving the corresponding linear system. Then the implicit surface can be reconstructed using the improved surface reconstruction method. In the process of surface reconstruction, all of the on-surface constraints are selected as the seed points. Figure 9 demonstrates comparison result of the voxel growing method with one and multiple seed points based on the on-surface constraints. The on-surface constraints are constructed by the sampling points. As mentioned, for the surface following method, the main problem is the incomplete reconstruction of sub-objects. Compared with the original voxel growing method, the improved method tends to recover the complete model as long as the sampling is sufficient. In the process of surface reconstruction, all of the on-surface constraints are selected as the seed points. Figure 9 demonstrates comparison result of the voxel growing method with one and multiple seed points based on the on-surface constraints. The on-surface constraints are constructed by the sampling points. As mentioned, for the surface following method, the main problem is the incomplete reconstruction of sub-objects. Compared with the original voxel growing method, the improved method tends to recover the complete model as long as the sampling is sufficient. In comparison with the marching cubes method, we test the performance of the improved voxel growing method and the leaf growing method on a variety of geological data sets, as shown in Figures 10 and 11. The surface following method does not modify the triangulation of voxels. Though the evaluation process of the three methods is different, the reconstruction results are the same when using the same size of reconstruction cell.  In comparison with the marching cubes method, we test the performance of the improved voxel growing method and the leaf growing method on a variety of geological data sets, as shown in Figures 10 and 11. The surface following method does not modify the triangulation of voxels. Though the evaluation process of the three methods is different, the reconstruction results are the same when using the same size of reconstruction cell. In comparison with the marching cubes method, we test the performance of the improved voxel growing method and the leaf growing method on a variety of geological data sets, as shown in Figures 10 and 11. The surface following method does not modify the triangulation of voxels. Though the evaluation process of the three methods is different, the reconstruction results are the same when using the same size of reconstruction cell.  In comparison with the marching cubes method, we test the performance of the improved voxel growing method and the leaf growing method on a variety of geological data sets, as shown in Figures 10 and 11. The surface following method does not modify the triangulation of voxels. Though the evaluation process of the three methods is different, the reconstruction results are the same when using the same size of reconstruction cell.  The performance of the reconstruction method mainly depends upon the number of interpolation constraints, the expansion order and the size of the reconstruction cell. The number of interpolation constraints determines the number of source points in the FMM evaluation. The expansion order determines the precision of evaluation, and the size of the reconstruction cell determines the model precision. Table 1 reports the timings of the evaluation of the marching cubes method, the voxel growing method and the leaf growing method on experimented examples. In the process of space division, the minimum bounding boxes of the sampling points are scaled by 0.01 times, respectively. Overall, because of the reduction of evaluation points, the voxel growing method has a better efficiency. The leaf growing method further improves the reconstruction efficiency. When the size of the reconstruction cell is large, the improvement in efficiency is not obvious. However, as the size of the reconstruction cell decreases, the performance advantage of the leaf growing method is more obvious. A series of experimental results show that the reduction of evaluation points can balance out the performance degradation for single point evaluation.

Conclusions and Discussion
In this paper, we improve an efficient implicit surface reconstruction method for the radial basis functions interpolant. The growing strategy of the surface following method combines both the evaluation and reconstruction processes. It ensures a smaller number of evaluation points and higher evaluation efficiency in surface following. The single point FMM is improved to evaluate points one by one. Based on the improved voxel growing method without manually specifying the seed points, the leaf growing method is developed to avoid a mass of redundant computation. The black-box FMM is utilized to evaluate the RBF interpolant with different kernels in a wide range of applications. The numerical results showed that one point evaluation in surface following obtains an optimum efficiency to the original approach. Compared with the existing methods, the improved method performs a better time and space efficiency.
The main limitation of the method is that the evaluations of the target points are no longer mutually independent because of the surface following process. It increases the difficulty of the parallel execution of the algorithm to a certain extent. To overcome this limitation, the space should be divided into several blocks. The reconstruction process between blocks does not affect each other. If the blocks are divided according to the octree structure in FMM, the evaluation process of each block can be also performed independently. Therefore, the leaf growing method is used to reconstruct the local implicit surface in parallel for each block. Finally, the results of each block are spliced into a complete surface.
To further improve the reconstruction efficiency and mesh quality, there are still some extensions that need to be studied. For example, the cells used in voxel growing are the same size, and the operations of single point FMM can be optimized for parallel implementation. To reduce the evaluation points, the function values and the gradient values can be combined in leaf growing to reconstruct the implicit surface adaptively.