Geological surface reconstruction from 3D point clouds

The numerical simulation of phenomena such as subsurface fluid flow or rock deformations are based on geological models, where volumes are typically defined through stratigraphic surfaces and faults, which constitute the geometric constraints, and then discretized into blocks to which relevant petrophysical or stress-strain properties are assigned. This paper illustrates the process by which it is possible to reconstruct the triangulation of 3D geological surfaces assigned as point clouds. These geological surfaces can then be used in codes dedicated to volume discretization to generate models of underground rocks. The method comprises the following: - Characterization of the best fitting plane and identification of the concave hull of the point cloud which is projected on it - Triangulation of the point cloud on the plane, constrained to the Planar Straight Line Graph constituted by the concave hull The algorithm, implemented in C++, depends exclusively on two parameters (nDig, maxCut) which allow one to easily evaluate the optimal refinement level of the hull on a case by case basis.


Method details
Historically, the process of generating numerical grids for geological applications is constrained to stratigraphic and fault surfaces. In fact, they are the main spatial elements that guide the zoning process of the model as well as its numerical discretization.
The paper describes a method developed to reconstruct these surfaces, represented in the space as generic point coordinates. Classes and functions are written in C ++ using the C library Triangle [3 , 4,10] , which operates on the plane. As an output it can provide a Delaunay triangulation constrained to a Planar Straight Line Graph , which, in this case, is defined by the concave hull of the point cloud. The Eigen library [5] , dedicated to linear algebra, was used for linear algebra operations.
The main steps of the methodology are schematized in Fig. 1 and described in detail in the sections which follow.

Best fitting plane projection
As already mentioned, the surfaces to be elaborated can be classified in stratigraphic and fault surfaces. If referring to a right-handed Cartesian coordinate system xyz where the z axis is vertically oriented, the stratigraphic surfaces are mainly orthogonal to the z direction, while faults are typically sub-vertical but do not have an orientation in space that can be determined a priori . As an example, the clouds of a stratigraphic surface (green) and of a fault surface (red) are shown in Fig. 2 .
Therefore, the xy-plane is a potential approximation plane for a stratigraphic surface, whereas, for a fault surface, it is necessary to calculate the principal components of the matrix consisting of the coordinates of the points of the cloud ( X ), determine the plane identified by them and project the points on it. The detailed procedure is the following: the point coordinates are expressed respect to a Cartesian reference system with origin in the barycenter of the cloud itself and the corresponding matrix X ∈ R 3 × N points is constructed. The eigenvalue matrix (D ∈ R 3 × 3 ) and the corresponding eigenvector matrix ( V ∈ R 3 × 3 ) of X X T were calculated using the Eigen library, then the projections were made on the plane identified by the first two components of the orthonormal basis obtained from V and corresponding to the two larger eigenvalues in D . The related functions are shown in Box 1 . Fig. 3 shows, by way of example, a cloud of points in space ( Fig. 3 A), the calculated BFP ( Fig. 3 B) and the projection of the points on it ( Fig. 3 C).

Hull "digging"
Once the reconstruction problem of the surfaces was brought back to a triangulation on the plane (BFP), the potential of the Triangle library was exploited, which implements a stable and efficient Delaunay triangulation algorithm and numerous refinement options. The syntax used for triangulation, without the imposition of any constraints, is shown in Box 2 , where the triangleIn and triangleOut structures are the input and output of the routine, respectively. Both are struct of the triangulateio type (see Supplementary Material section for details), whose declaration and definition can be found in the files triangle.h and triangle.c of the Triangle library.
However, if not constrained to the boundary nodes, the point cloud triangulation on the convex hull might lead to a poor quality discretization when projected back in the space. This is the reason why a series of triangles (red area of Fig. 4 ) can be observed, which are incorrectly defined because  they connect points that should be tagged as boundary nodes. As a consequence, the triangulation of the fault surfaces in space shows a series of erroneous "folds". Therefore, the need to more accurately identify the polygon that describes the boundary of the BFP projection of the point cloud, i.e. the concave hull , became evident. Routines based on the concept of αshapes, available in the CGAL library [6 , 7] , were preliminary tested, but unsatisfactory results

Box 1
Functions for deriving the BFP e consequent projection. The points are saved in the Eigen::MatrixXd xyz and their projection on BFP in xyProj.

Box 2
Triangle syntax for triangulation of the 2D point cloud without boundary constraints (see details on Triangle switches in Supplementary Material section) .    5. Sketch of the exploration step of the "digging" algorithm: nodes classified as internal in white, nodes and segments belonging to the boundary polygon in red. In a hypothetical iteration, the p 1 p 2 segment is tested. In particular, the p 1neighborhood is defined by the circumference with radius r and center p 1 and the algorithm searches for internal nodes that fall within it.
were obtained due to the different concavities characterizing the point cloud boundary. For this reason, taking as a reference the "dig criterion" of the convex hull introduced by Park and Oh [1] , we implemented a variant of the algorithm that is able to optimize the exploration of both the extension and the shape of the neighborhood (for the search of possible boundary points). The algorithm is based on the idea of incrementally "digging" the convex hull by iteratively adding a node previously classified as internal to the boundary polygon. If the node passes the admissibility test, it is inserted into the appropriate position and the existing segment of the polygon is "broken". Fig. 5 shows the exploration scheme of the "digging" algorithm. The starting polygon of the algorithm is the convex hull obtained as the output of Triangle and the depth to which the hull is "excavated" is determined by the parameter indicated with nDig ∈ [0,1]: it defines the width of the neighborhood exploration through the relation: r = edgeLength * dig. The neighborhood to explore is wider when the dig parameter approaches 1. In a hypothetical iteration of the algorithm, the edge p 1 p 2 of the current boundary Box 3 hull::updateHullDig(double dig) -pseudo-code of the function for the construction of the concave hull using the "digging" strategy. polygon is tested. In particular, the p 1 -neighborhood is defined by the circumference with radius r and center p 1 and the algorithm searches for internal nodes that fall within it. In Fig. 5 , p neigh (yellow) is selected and it is tested, i.e. it is verified that the trial boundary edges p 1 p neigh and p neigh p 2 do not intersect any edge of the current boundary polygon and that no internal points fall in the triangle p 1 p neigh p 2 . If the point p neigh passes the test, then it is tagged as boundary node and the p 1 p 2 edge is "broken". The Box 3 reports the pseudo-code instructions of the updateHullDig function and Table 1 summarizes the main variables and type and provides their description. The additional routines which are called are reported in Box 7-Box 17 (see the Supplementary Material section ) for completeness. Fig. 6 shows an example of evolution of the "digging" algorithm. The initial convex hull (in red) and four concave hulls of increasing number of nodes (from orange to yellow) are represented. The result is the green boundary.
In the development of the method, some distinctions had to be made between fault surfaces and stratigraphic surfaces due to the intrinsic characteristics of the latter. First, we observed that the calculation of the BFP was not needed as, generally, a BFP can be approximated by a xy-plane of a right-handed Cartesian coordinate system with z as the vertical axis. Secondly, the point clouds that we have to process are typically surfaces deriving from geological modeling and sampled with a constant step on the plane. It follows that the points, once projected onto the plane, have a regular distribution and, in particular, they are aligned with the coordinated axes or with the diagonals, falling within the limit configurations that make the implemented crossing algorithm less accurate [8] . Indeed, one of the implemented core function belongs to the family of strategies called point-inpolygon , whose criticalities are extensively tested in [9] . Table 1 Description of the main variables defined in the pseudo-code.

Table 2
Reconstruction of fault surfaces with the updateHullDig routine. The number of the points in the cloud, the values assigned to the dig parameter, the characteristics of the triangulation referred to the initial convex hull and the values relative to the final concave hull are given.      In order to overcome these limitations, we decided to exploit the information deriving from the preliminary triangulation, i.e. the triangulation delimited by the convex hull , so as to apply a different "digging" strategy. We observed that, since the sampling is regular, the quality of the mesh, measured as edge ratio (i.e. ratio between the maximum edge length and the minimum edge length of each triangle of the grid), is roughly constant over the whole domain. Triangles characterized by edge ratio higher than average are found exclusively close to the boundary polygon where it is known that the  convex hull is not adequate to describe the profile of the point cloud. As an example, Fig. 7 shows a portion of the triangulation of the point cloud of a stratigraphic surface projected onto the xy-plane, delimited by the convex hull (nodes in red). The color map describes the trend of the edge ratio. Triangles that deviate from the average value (approximately equal to 1.2) are just the triangles that should be excluded from the triangulation because they are external to the polygon that describes the boundary of the projected cloud (nodes in green). Based on such considerations we decided to calculate the mean ( t . mu ) and the standard deviation ( t . sigma ) of the triangle edge length distribution of the initial discretization. At the same time, two lists of triangles were populated: the first, named outlierTriangles, contains all the triangles that have at least one edge of length greater than thresholdLength = t. mu + maxCut * t . sigma , where maxCut ∈ R + is a parameter set by the user; the second, named starterTriangles, is a subset of the previous one and includes the triangles that have one edge in common with the starting convex hull . Therefore, the second list consists of those triangles that should be potentially excluded from the triangulation to improve the definition of the hull. Once the two lists are initialized, the algorithm (whose steps are detailed in Box 4 ) iterates over the elements of the starterTriangles. In the representation of Fig. 8 A, the list consists of two triangles: T 1 and T 2 . In particular, both edges s 1 e 1 and s 2 e 2 belong to the current boundary polygon. Thus opposite vertices p i ( i = 1,2) are tested in order to be tagged as boundary nodes and to be inserted in the concave hull , i.e. the algorithm verifies that the trial boundary edges s i p i and p i e i do not intersect any edge of the current boundary Box 4 hull::updateHullWipe(triInfo t, double maxCut) -pseudo-code of the function for the construction of the concave hull on the plane. The strategy eliminates triangles by evaluating the length of their edges and whether they belong to the boundary polygon.

Box 6
Function for the transformation to the original reference system from local coordinates on the plane. polygon and that no internal points fall in the triangle s i p i e i . If the admissibility test is passed, T 1 and T 2 are removed and the boundary polygon updated ( Fig. 8 B). Eventually, we go through the triangles that have a side in common with the removed triangles exploiting the neighborlist (see Table 1 ). If they belong to the outlierTriangles list, they too are added to the starterTriangles list. In the example in Fig. 8 B, triangle T 3 is added.
The algorithm continues to update the starterTriangles list until it is empty. The resulting boundary polygon will then be the final concave hull .
The methodology involves a second triangulation on the plane where the Planar Straight Line Graph defined by the identified polygon is imposed as a constraint.
The corresponding string used to run the Triangle library is given in Box 5 . At last, the constrained triangulation is re-projected in the 3D space to obtain the reconstructed surface ( Box 6 ).
The resulting triangulated surface inherits the IO structure of the Triangle library, i.e. the triangulateio C-struct, whose pointlist array is updated with the 3D re-projected coordinates. It is observed that during the triangulation process no points are added or deleted thus there is a perfect correspondence between the input and output points of the cloud.

Method validation
In order to illustrate the effectiveness of the presented method, the reconstruction processes of 2 fault surfaces and 2 stratigraphic surfaces are reported below. Table 2 shows the details of the application of the updateHullDig routine to the cases named Fault 12 ( Fig. 9 , Fig. 10 ) and H53 ( Fig. 11 ,    Table 3 we refer to the cases called Surf Top ( Fig. 13 , Fig. 14 and Fig. 15 ) and Erosional 50 ( Fig. 16 , Fig. 17 and Fig. 18 ) where the updateHullWipe routine was applied.
We observe that in all the presented cases the triangles of the initial triangulation (constrained to the convex hull ) that generated erroneous "folds" (in red) are correctly removed from the final triangulation (constrained to the concave hull ). Removed triangles originally connected points that are classified as boundary nodes in the final hull.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.