The control of rock slope movements in open-pit mines is important for the secure work and the continuity of production (Nunoo et al. 2016). Any uncontrolled rockslide or landslide may have several impacts on the safety and economics of an open-pit operation. Thus, the monitoring of pit walls at open-pit mines is a crucial task to manage the stability of them (Nghia and Maciaszek 2010).

There are several slope monitoring methods divided into three categories: (1) visual inspection, (2) surface measurement and (3) subsurface measurement (Nunoo et al. 2016). While for visual inspection, mine operations personnel working in the pit watch for any potential unsafe pit wall behavior (Nunoo et al. 2016), surface and subsurface measurements use complex technologies or special instrumentation in performance. For surface measurement, surveying devices like global positioning system (GPS) receivers (Lipecki 1999; Younger 2001), total stations and prisms (Lazzarini 1977) are widely employed with high accuracy. However, these techniques gathers points of objects only, therefore, these are time consuming and not cost-effective methods (Abellán et al. 2009).

Recent advances in the field of geospatial technologies has provided several new monitoring techniques for the detection of rock slope instabilities across a wide range of spatial and temporal scales such as Interferometry Synthetic Aperture Radar (InSAR) (Amelung et al. 1999; Krawczyk et al. 2012), space-born Light Detection and Ranging (LIDAR), and terrestrial LIDAR. All three techniques can measure whole visible surface of objects, but the first two techniques allow only to collect data for fixed periods when satellites or planes has flown above, whereas the terrestrial LIDAR is a rapid method and it can be used to measure large areas (Bitelli et al. 2004; Ferraz et al. 2016; Tang et al. 2014). In addition, terrestrial LIDAR can replace tacheometric measurements in different tasks, such as: structures’ health monitoring (Park et al. 2007), rivers (Williams et al. 2014) and forest measurements (Simonse et al. 2003), and creating geological documentation (Buckley et al. 2008). Therefore, the terrestrial LIDAR has been proposed for topographic mapping in open-pit mines (Barbarella and Fiani 2013; Jaboyedoff et al. 2012). The main advantage of the terrestrial LIDAR is that the accuracy of the derived products is very high (Boehler et al. 2003). In addition, temporal resolution can be decided by surveyors to measure all visible surface of the object. However, the accuracy assessment of topographic surfaces in open-pit mines derived from the terrestrial LIDAR is still rare. Moreover, the terrestrial LIDAR has not been studied yet as extensively as other methods like GPS or total station measurements.

There are several different methods for converting measured point clouds into a 3D polygonal (Cignoni et al. 1998; Remondino 2003b), including methods for generating irregular network, rectangular grid and the contour overlap (Yao et al. 2014). The reconstruction of precise surfaces from unorganized point clouds derived from laser scanner data is a very hard problem(Remondino 2003b). In addition, there is a lack of comparative studies of different interpolation algorithms for generating topographic surface at open-pit coal mines. Surface reconstruction by mesh algorithms allows to create continuous representation of the object based on discrete data. In other words, the information about object structure is known (Bucksch and Lindenbergh 2008). Open-pit mines’ areas are usually full of bed rocks, flat stones and sharp edges. Owing to that, it is easy and accurate to present them as a mesh surface. Nevertheless, using a non-uniform rational basis spline (NURBS) surface pretend to be appropriate in this case due to the possibility of the assumption of the pit mine area as a continuous surface.

This work addressed this issue by comparing three different algorithms based on the Delaunay triangulation, namely Delaunay 2.5D XY Plane, Delaunay 2.5D Best Fitting Plane, and Mesh from Points for the reconstruction of surfaces from terrestrial LIDAAR data. It also presents the application of the polygonal mesh to determine displacements. The mesh surface was generated from the point cloud data. It can be used as a reference object for further measurement series as well as a basis for generation of a NURBS surface. Moreover, quality of created surfaces, mesh quality improvement approaches and possible sources of errors were discussed in article. Additionally, the application of NURBS surfaces generated from mesh networks in modelling the noisy point clouds was present, and their accuracy were assessed based on two independent point clouds.

The point clouds used for the research were obtained at the Akademia Górniczo-Hutnicza (AGH) experimental mine and at the Coc Sau open-pit mine. The AGH experimental mine is a museum and a laboratory placed at the AGH University of Science and Technology (AGH UST-https://www.agh.edu.pl/en) in Cracow (Poland). It is unique research and teaching place, presenting the complete infrastructure needed in the underground coal exploitation. For this area, 3D point cloud was obtained by using Leica C10 (Leica Geosystems 2011). The main characteristic of this terrestrial LIDAR instrument was summarized in Table 1. The obtained point cloud after being merged and cleaned by Leica Cyclone software was present in Fig. 1. The density of the point cloud was around 2500 points per 0.01 m² and the points were distributed uniformly.

Point cloud of the test field wall from AGH experimental mine

Coc Sau open-pit mine, located in the Mong Duong Ward, Cam Pha city, Quang Ninh province, Vietnam, around 200 km to the east of Hanoi city, was selected as second case study. Recently, the mine has been exploited at a depth of − 200 m above sea level. The mine is using modern technologies including drilling explosion, coal excavation, transportation, mineral processing and sorting to consumption (Khan 2017). For this area, 3D point cloud for the topographic surface was scanned using GeoMax Zoom300 LIDAR (Fig. 2), the main characteristic of this terrestrial LIDAR instrument was summarized in Table 1.

A photo of the Coc Sau open-pit coal mine (the photo was taken by Michał M. Buczek on November 7, 2016)

In Coc Sau case study, the data was collected from one station in local coordinate system. Owing to that, there was no need to use control points. The point cloud was processed and cleaned in XPAD program. According to the user guide, the accuracy of the point cloud should be around 6 mm per 50 m (GeoMax AG 2015). The density of the obtained point cloud was around 20 points per 0.01 m². From the post-processed point cloud, a small part presenting the highwall slope was chosen (Fig. 3). Due to mine conditions, the point cloud was obtained from far stations, slopes were scanned at a sharp angle and the obtained data was noisy; therefore, this dataset was chosen to assess the practical application of the tested models on the real-life with low quality data. This part of the data set was used for all analyses and was called as “the original point cloud” (OPC).

The point cloud of the chosen highwall from Coc Sau open-pit mine

3.1 Data preparation

Before mesh generation, pre-processing operation must be performed. Firstly, the data set should be sampled. After that noise reduction, outliers excluding, and holes filling should be carried out to remove errors. Next step of creating polygonal mesh surface is global topology determination. Then, the polygonal surface can be created. Post-processing is a part of the model construction, when an object surface is editing automatically or manually—corrections are performed, mesh could be repaired, smoothed or split into parts (Remondino 2003b).

3.2 Mesh generation

3.2.1 Background of the mesh algorithms

Mesh generation was based on Delaunay triangulation, which was found in 1934 as a dual graph of Voronoi diagram (Delaunay 1934). It is a common tool used in graphic and 3D modeling applications to convert point cloud into 3D polygonal model (Fig. 4).

Voronoi diagram (left) and Delaunay triangulation (right) of the same set of points (Remondino 2003a)

2.5D Triangulation is the process based on Delaunay algorithm, which is used as a interpolation function. For set of points P and real, unique elevation function f(x,y) at each point p(x,y) ∈ P the linear interpolating function F is created. F(p) is weighted average of the elevation of triangle’s vertices which contains p (Remondino 2003a).

In Micro Station V8i triangulation method uses horizontal plane for data computation and determination of connections between network nodes.

Delaunay 2.5D XY Plane and Delaunay 2.5D Best Fitting Plane are algorithms available in CloudCompare. First algorithm is based on the point cloud projection on the 2D XY plane. After projection 2D points are triangulated and new-build mesh structure is used to create 3D surface (Girardeau-Montau 2018a). The difference of the second one is that point cloud is projected on the best fitting plane not to XY plane. In this process, least square method is used (Girardeau-Montau 2018b).

3.2.2 Case study in AGH experimental mine

In the AGH experimental mine a part of the longwall face was chosen. The wall’s surface is almost a planar with shallow pockets. A part of the point cloud with visible pockets was illustrated in Fig. 1. The distance between LIDAR station and the longwall was less than 4 m. The obtained pointset was chosen to evaluate the accuracy of fitted mesh networks and NURBS surfaces. This high-quality point cloud was used to determine the accuracy of the surface models. The meshes were generated using three different algorithms available in MicroStation software and CloudCompare software: (1) MicroStation’s Mesh from Points—Global Z Direction—(AGH MS); (2) CloudCompare’s Delaunay 2.5D (XY plane)—(AGH CC); (3) CloudCompare’s Delaunay 2.5D (Best Fitting Plane)—(AGH BF). These meshes were called ‘original mesh’. After that, the results were rectified and optimized, and called ‘optimized mesh’. The worklow for the AGH experimental mine case study was illustrated in Fig. 5.

The AGH experimental mine case study workflow

3.2.3 Case study in Coc Sau open-pit coal mine

Similar to AGH case study, the mesh networks were generated from the OPC using different algorithms available in MicroStation software and CloudCompare software: (1) MicroStation’s Mesh from Points—Global Z Direction—(MS mesh); (2) CloudCompare’s Delaunay 2.5D (XY plane)—(CC mesh); (3) CloudCompare’s Delaunay 2.5D (Best Fitting Plane)—(BF mesh). However, as mentioned in Sect. 2, the OPC of Coc Sau was quite noisy, which could influence the quality of meshes. Therefore, the comparison with data not used to create networks could be more accurate assessment. Additionally, in order to assess the practical application of the tested models on the real-life with low quality data, the second point cloud was generated to simulate the process of displacements. To do this a disturbance smaller than 10 cm was added into the OPC. Using the propagation of uncertainty, authors computed the maximal displacement of single point in 3D around up to 173 mm (approximate to 100 \(\sqrt 3\) mm). The matrix of disturbances was generated with values from the range between − 0.1 m and 0.1 m with normal distribution and the created point cloud is called “the disturbed point cloud” (DPC). The worklow for the Coc Sau open-pit mine case study was illustrated in Fig. 6.

The Coc Sau open-pit mine case study workflow

3.3 NURBS surface generation

Generated meshes may have errors such as intersections, wrong normal, flipped triangles. Therefore, the mesh rectification and optimization are necessary before any further usage. Then, the rectified mesh was used to generated NURBS surfaces (Krishnamurthy and Levoy 1996), which is a mathematical model for generating and representing curves and surfaces (Domingo et al. 1995).

The NURBS surface approximates the input mesh network and a decrease of accuracy is inevitable. On the other hand, the shape of the mesh network depends on the pointset and their errors. An example is present in Fig. 7. The pointset (blue color) with local irregularities was modeled with a mesh network and a NURBS surface. For both models, the disturbed place influences the models’ plot. The deformation was not as big as it was in the pointset. If the pointset has random irregularities, the replacing it with models may reduce its noise.

Irregular part of point cloud (blue) modeled with the mesh network (red) and the NURBS surface patches (black) from two perspectives

The shape of NURBS surface is determined by its order, a set of weighted control points, and a knot vector. It can be built from numerous patches (black lines in the Fig. 7). The surface patch is defined by the mathematical formulas: interpolated (1) and approximated (2):

where:\(t = \left\{ {t_{0} , \ldots , t_{n} } \right\},u = \left\{ {u_{0} , \ldots , u_{r} } \right\}\): function knots, \(d_{i,j}\): control points, \(N_{i} ,N_{j}\): basis functions, n, r: number of control pointsComputer program GeoMagic Design X was used to compute the NURBS surface from optimized meshes.

3.4 Accuracy assessment

To assess the accuracy of the results, three parameters were used. The first one, the maximal distance, is the distance between the furthest point from the point cloud and the surface. The second one, the average distance, is the mean distance between all points from the point cloud and the surface. The last one is the standard deviation.

The accuracy was assessed from the distance statistics for combinations of point cloud—mesh and mesh—mesh. The statistics were calculated for the pairs as follows: (1) original mesh network—original point cloud; (2) optimized mesh network—original point cloud; (3) original mesh network—optimized mesh network; (4) original mesh network—disturbed point cloud; (5) optimized mesh network—disturbed point cloud.

The same statistics were used to assess the NURBS surface. The distances between the point sets and the NURBS surfaces were computed. The point cloud was also divided into few ranges according to distance to the surface and distance maps were created.

4.1 Mesh networks

By using the method mentioned in Sect. 3.2 the meshes for AGH and Coc Sau study case were generated. The accuracy of fitting all mesh networks into point cloud was present in Tables 2 and 3. CloudCompare software could have a problem with generated original meshes, as a result all the statistical parameters were equally 0.

For each optimized mesh network in AGH case study, the average distance was lower than 2 mm, and the deviation lower than 7 mm. Those results are similar or even smaller than the accuracy of the single measurement of the Leica C10 laser scanner, which are 4 mm for position and 6 mm for distance (Leica Geosystems 2011). Thus, the mesh optimization process has preserved the accuracy within the accuracy of the point cloud.

For Coc Sau open-pit mine, the results showed that Mesh from Points algorithm created higher quality mesh network than the other ones. Rectifying tools from Autodesk Netfabb (Repair) and GeoMagic Design X (Mesh Doctor) have not found any errors, holes, and discontinuities inside the MS mesh, in contrast to the CC mesh and BF mesh. The Best Fitting Plane algorithm created networks with 10 times more triangles than the XY plane algorithm. It was caused by the irregular distribution of the points in the point (Fig. 8).

A part of highwall at Coc Sau modeled by two algorithms: CC mesh (left) and BF mesh (right)

Analyzing the BF mesh, it can be notice that its structure depends on the distribution of points in the point cloud, and the CC mesh better present the real shape of the scanned slope. However, both mesh networks required a rectification and all the errors were removed in the iteration process.

In the next step, the mesh networks were optimized. As a result, the new meshes were created with regular triangles of the similar size. This process also lowered the density of triangles about 3–4 times. Figure 9 illustrated the same part of the BF mesh before and after the optimization. As can be seen, the surface was smooth, and the characteristic lines were still visible. The accuracy of fitting all mesh networks into point cloud were present in the Table 3. There was no similar correlation to the one for the AGH test site present in Table 2. However, it is worth notice that after optimization CC mesh and BF mesh were quite similar. The average distance for optimized MS mesh was 4 times bigger comparing to that of optimized BF mesh.

Comparison of the chosen part of the mesh CS-BF before (left) and after (right) optimization process

4.2 Accuracy assessment of the optimized meshes

To assess the accuracy of the optimized meshes in comparison to the original meshes, the same three parameters were used. Table 4 showed the results of meshes comparison before and after optimization. The MS mesh did not have any errors to repair, but still there was a need to recreate triangles to have similar size and similar density in the whole network. The average distance for the CloudCompare’s methods was under 1 mm, while for the MicroStation it was dozens of time bigger. Analyzing deviation values, the original and optimized CC mesh can be assumed to be the best fitted mesh.

As can be seen in Table 4, the average distances between original and optimized meshes are equal to 0.1 mm ± 2.9 mm for CC mesh, and to 0.3 mm ± 8.2 mm for BF mesh. The process of optimizing mesh network decreased the precision of fitting surface into points. The precision loss might depend on the shape of input mesh network, as well as the dispersion of the input pointset. Additionally, all meshes had great values of maximal distance. This could be caused by holes and single outlier points. Figure 10 showed the distance map for the BF mesh before and after optimization. It showed that the triangles generated for the holes indicated by red color had bigger distances.

Deviation between original mesh BF and optimized mesh BF

4.3 Accuracy assessment based on DPC

Summary of the accuracy assessment for the new point cloud, which was not correlated to generated meshes, was shown in Table 5. Once again, all parameters for MS mesh were few times bigger than the other meshes. For CC mesh and BF mesh, average distance was lower than 1 cm, and after optimizing it was lower than 2 cm. However, the deviation was about 3–4 times bigger than the average distance.

It can be seen that the maximum distances were greater than the noise value (173 mm) added to the OPC. The average distance for all tested objects was smaller than 43.3 mm, that is ¼ of the maximum displacement error. The standard deviation in all cases was lower than 86.6 mm, that is ½ of the maximum error. Thus, the precision was on the level of 2P = 95%.

4.4 NURBS surfaces

The optimized mesh networks were used to generate NURBS surfaces. Figure 11 depicts the part of the AGH MS mesh and the corresponding part of the NURBS surface. The accuracy of the surface was shown in Table 6. It can be seen that the average distance was 4–5 times bigger than that of the mesh networks. The values were closed to the laser scanner accuracy. Thus, the final surface may be used as a model of the point cloud on the known level of accuracy.

A part of the AGH MS mesh network (left) and the NURBS surface patch (right)

In this case, the test field is quite flat. Therefore, the model of a planar surface was fitted into the pointset. The result showed that the average distance was about 5 times bigger than NURBS models. It might be caused by the noise as well as the shallow pockets on the wall. It proves that simple geometric objects cannot be used to model a noisy point cloud with high accuracy, even in simple cases.

The same approach was applied on the optimized mesh networks of Coc Sau mine. The rectangular grid with size of 2 m× 2.5 m was generated from the spline curves. The spline curves could not have discontinuities and loops. The joined splines created the wireframe build of quadrangles. Each rectangular of the wireframe was used to compute one patch of the NURBS surface using Legendary Boundary Fit tool. The grid network (wireframe) and the NURBS surface was shown in Fig. 12.

NURBS surface fitted into optimized mesh CC

The average distance was computed for each NURBS surface, and the results were presented in Table 7. Same as in AGH test field, the values were 4–5 times bigger than accuracy of the optimized mesh networks. The quality of MS mesh was the worse, what strongly influenced the NURBS accuracy. The point cloud cannot be assumed as a flat surface, therefore the free-form surface was computed using RhingResurf’s Single Surface from Points. The accuracy was worse than for the NURBS surfaces computed with mesh networks. Therefore, the presented approach gave better results than approach modeling the raw point cloud. Additionally, the results proved that loss of accuracy was within the accuracy of the input point cloud.

The NURBS surfaces were generated from the optimized meshes created from the OPC. For each surface, the distance map for the DPC was created. The example visualization for BF mesh was shown in Fig. 13. Five groups of distances were classified: (1) distance smaller than 43 mm; (2) distance from 44 mm to 86 mm; (3) distance from 87 mm to 131 mm; (4) distance from 132 to 173 mm; (5) distance greater than 173 mm. Statistic of the distances by groups were present in Tables 8, 9 and 10.

BF NURBS surface covered by disturbed point cloud (deviation between data sets is shown by colors: blue—group (1), light blue—group (2), green—group (3), orange—group (4), red—group (5)

As shown in Tables 8, 9 and 10, the best NURBS surface was generated from CC mesh. Almost 50% of all points had distance smaller than ¼ of the maximum distance error, and almost 75% of points had distance ½ of the maximum error. Only single points were further than maximum distance of disturbance. The result for the surface from BF mesh was slightly worse.

On the other hand, the worst results were for the surface from MS mesh. Only one-third of the points had distance smaller than ¼ of the maximum error. There was a big group of points (almost 6%) laying further than maximum error. The MicroStation generated the best raw mesh network, without any errors. However, the quality of fitting the original and the disturbed point set was the worst.

The processing of a raw point cloud is a resource-intensive process. In addition, it cannot be applied directly in CAD designing. Thus, one has to convert point cloud into simpler 3D graphics objects like polygon meshes and surfaces. The presented methods allowed creating digital models of the real world object from the laser scanning data. The networks optimization allowed creating smooth NURBS surfaces. These tools can be especially applied with the noisy datasets instead of simple geometrical objects.

The NURBS surface is continuous and smooth, and allows generating a map of distances for all elements from the point cloud. The laser scanning gives information about all measured object, not only the chosen control points. Thus, the presented approach may be an alternative to the conventional methods, giving more complete and accurate data about observed object.

Applying the NURBS surface to approximate the point cloud reduces the influence of noise and single, dispersed points on the model’s shape. The NURBS fitting accuracy into pointset is 4–5 times worse than optimized mesh. It is the same for AGH and Coc Sau point clouds. Therefore, the quality of the point cloud has not influenced the quality of NURBS surface. It depends strongly on the quality of the input mesh network.

Further work on this topic will focus on improving the quality of the laser scans, to improve the meshes. They also want to replace the second series with real observations of displacement of the unknown value. The comparison with classic methods would give the final answer if this method is enough precise to determine deformations.