STOCHASTIC SURFACE MESH RECONSTRUCTION

A generic and practical methodology is presented for 3D surface mesh reconstruction from the terrestrial laser scanner (TLS) derived point clouds. It has two main steps. The first step deals with developing an anisotropic point error model, which is capable of computing the theoretical precisions of 3D coordinates of each individual point in the point cloud. The magnitude and direction of the errors are represented in the form of error ellipsoids. The following second step is focused on the stochastic surface mesh reconstruction. It exploits the previously determined error ellipsoids by computing a point-wise quality measure, which takes into account the semi-diagonal axis length of the error ellipsoid. The points only with the least errors are used in the surface triangulation. The remaining ones are automatically discarded.


INTRODUCTION
Terrestrial laser scanners (TLSs) capture the geometry of target object or scene in the form of dense point clouds.Any point in the scan data is contaminated by the random errors.These errors propagate through the steps of data processing, namely preprocessing, co-registration, mesh generation and 3D model reconstruction.Estimating the random error pattern of every individual point is essentially important for 3D model related tasks.The range of applications are various such as surface matching (Akca and Gruen, 2005;Gruen and Akca, 2005;Akca and Gruen, 2007), 3D object modelling and surface mesh generation (Akca et al. 2006a(Akca et al. , 2006b;;Akca et al. 2007), and surface comparison (Zhang et al. 2006;Baltsavias et al. 2007).
The target object has to be scanned from multiple standpoints to ensure full coverage.Such a scanning configuration will result in overlaps between the consecutive scans.In this case, overlapping parts of the target object will be sampled by two or more point clouds each of which have different error patterns.The most common solution of this redundancy problem is the subsampling the data and then generating the 3D surface mesh.This solution is not optimal, since the random error characteristic of each individual point is anisotropic.We aim at developing a volumetric resampling method in which the low quality points are discarded and the high quality points are used in the subsequent modelling steps.Such an approach will yield 3D surface models with the least errors.Contributions of this paper can be summarized as follows: • Development of a practical and generic point error model for TLS-derived data • Utilizing of a stochastic surface mesh reconstruction method based on the developed point error model.

Related work on (general) mesh reconstruction
Significant amount of works have been done on surface mesh reconstruction.These works can be classified to four main groups (Santos et al., 2012): Delaunay-based, surface-based, parametric surfaces-based and volumetric methods.
Delaunay-based methods generally use unorganized point clouds as the input, allowing a Delaunay complex to be established.Using the Delaunay complex, mesh reconstruction can be realized by means of alpha-shapes (Bernardini et al., 1999a).An often used method Crust reconstruction was presented in (Amenta et al., 1998) which is based on Voronoi diagram.This method was further improved using Cocone based reconstruction methods (Dey and Goswami, 2003).
Surface-based methods take the surfaces of each scan or partial surfaces of the target object as the input.As a result of integration, a single triangular mesh (or surface) is reconstructed.One of the most known methods was presented by Turk and Levoy (1994) which zippers the meshes in order to obtain a consensus surface.The ball-pivoting algorithm (Bernardini et al., 1999b) can also be considered as a surfacebased reconstruction method.It incrementally interpolates the surface triangulation.
The work of Sclaroff and Pentland (1991) can be considered as pioneering in the parametric surface-based works.They developed an approach which is based on the finite element model (FEM).Kazhdan et al. (2006) considered reconstruction as a spatial Poisson problem.This approach is improved by the implementation of a parallel computing algorithm (Bolitho et al., 2009).For the solution of the Poisson problem, a screening term is introduced in addition to the mathematical framework.So far this effort increases the quality of the integration (Kazhdan and Hoppe, 2013).
The volumetric methods implicitly use the voxels which is a kind of 3D cuboid data structure where each voxel has attribute values.Hoppe et al. (1992) developed an approach which first generates a signed distance function (SDF) from the unorganized points using a Euclidean minimal spanning tree, and then applies the marching cube algorithm (Lorensen and Cline, 1987) to reconstruct the surface mesh.Curless and Levoy (1996) developed the volumetric range image processing (VRIP) algorithm which initially creates a dense volumetric grid.For each voxel, the weighted signed distance of points to the nearest range surface is computed in a line of sight direction.

Related work on error consideration
The surface reconstruction methods relying on the moving least squares (MLS) approach have been used both suppressing the noise and reconstructing local surface elements.Such an algorithm presented by Jones et al. (2003) used the original surface mesh as input.A local weight function was introduced using the spatial neighbourhood of each vertex.Then, position of each vertex was altered with respect to the local weight function.This algorithm not only de-noises the surface but also preserves sharp details.Point clouds can also be used as input.
For example in Pauly (2003), local surface analyses were carried out with the MLS, followed by computation of SDF for surface mesh generation.
Just after the acquisition step, the point cloud data may be preprocessed to eliminate noise and other errors.The preprocessing step is generally required in case of low density sampling, holes and surface reflectance induced problems.Weyrich et al. (2004) presented a tool box implemented based on the MLS approach.It can be used to remove outliers, to fill holes and to smooth the point cloud data.Thereafter, the resulting point cloud can be used to generate high quality surface meshes by employing standard surface triangulation methods.
The study presented in Adamson and Alexa (2006) differs from the others by introducing an anisotropic MLS method.The anisotropy is based on the individual weight functions defined for each individual point, rather than a point and its spatial neighbourhood.In computation of weight function, only the principal curvatures are used.Each point is associated with its corresponding ellipsoid coming from covariance analysis of the weighting function.
In case of sparsely sampled point clouds, the local LMS methods may lead insufficient results.To overcome this issue higher order algebraic surfaces such as sphere can be fit (Guennebaud and Gross, 2007).This method generates successful surface meshes from sparsely sampled data containing high curvature details.The MLS concept can be used with the robust statistical methods for surface mesh reconstruction (Oztireli et al., 2009).This study has advanced the MLS by integrating the non-linear Kernel regression method.It does not require any pre-processing step.The sharpness of features can be controlled by the user.It has realtime reconstruction capabilities.
Such a large number of studies show the relevance of the problem.A fully satisfying solution, which mathematically formulates the physical nature of the instrumental and environmental errors and applies it to the surface modelling tasks, has to be still designed, realized and justified.The point cloud data has certain error distributions.Its pattern can be investigated by means of positional uncertainty which is influenced by multiple parameters.Angular (mechanical) stability, sensor-to-object distance, incidence angle of the incoming laser beam, and surface reflectivity are the most significant ones.It results in a heteroscedastic (point dependent), anisotropic and inhomogeneous point error distribution.This fact has to be considered cautiously in the surface mesh reconstruction.The contribution of each induvial point to the final surface has to be evaluated separately.We developed such a surface mesh generation method.The point error modelling is given in the second section.The subsequent surface generation method is given in the third section.The experimental results are presented in the fourth section.

ANISOTROPIC POINT ERROR MODEL
The TLS systems operate in a spherical coordinate system measuring the range (ρ), vertical (α) and horizontal (θ) angles as the direct observations.The Cartesian coordinates of any -th point are computed from the spherical observations = as: The spherical observations can be computed by reversing Equation (1) unless they are provided by the TLS software.Using the law of error propagation, the variance-covariance matrix ∑ of any point in the Cartesian coordinate system is: where is the Jacobian matrix of the partial derivatives of the Cartesian coordinates with respect to the range ( ), vertical ( ) and horizontal ( ) angle observations.∑ is the a priori variance-covariance matrix of any observation point where, , " and # represent a priori variances of the range ( ), vertical ( ) and horizontal ( ) angles, respectively.Nondiagonal elements of this matrix are zero under the assumptions that (a) the TLS is well calibrated so that there is no sensor caused distortion and (b) there is no physical correlation among the direct observations.
The parameters of the error ellipsoid can be calculated from the principal components of the variance-covariance matrix ∑ as shown in Equation (4).The error ellipsoid fictitiously represents the magnitude and direction of the random error of the associated point.At this stage, the error ellipsoid of every individual point can be computed provided that the ∑ in Equation ( 3) is known.In order to fill the ∑ matrix up, a priori precision values , " and # should be computed.We developed a practical methodology to compute them, which is explained in the following.
2.1 Computation of the angular precision values , -and , .
A static and repetitive measurement configuration is proposed.The vertical ( ) and horizontal ( ) angle precisions are determined by the repeated scans of the same environment when the TLS is set up firmly static.It is expected that each point should coincide with its conjugates in the other scans.The deviations of the conjugate points are relevant to the angular repeatability of the scanner system.
In the experimental studies, the same environment was scanned five times from the same station.The five conjugate points of the same laser ray were selected and corresponding vertical ( ) and horizontal ( ) angles were computed.The root mean square error (RMSE) values of their discrepancies were computed both for the vertical and horizontal angles.These RMSE computations were repeated at least four laser rays which towards the four main directions North, East, South and West or equivalents.The mathematical means of each of these four RMSE values yield the a priori angular precisions " and # , respectively.

Computation of the range precision value , /
We have developed an empirical formula for the range precision ( ) computation based on fieldwork and validation studies given in our previous work (Ozendi et al., 2017).It takes into account the distance between the TLS and the point, incidence angle of the incoming laser beam, and reflectivity of the object surface.
= 0 1 2 1 3 4 567 8 (5) and 9 : = ; < = > , for : B : C 0 , elsewhere where <, >, H and I are the coefficients, is the scanner-to-point distance, J is the incidence angle and : is the intensity value of the point.If provided by the associated software of the scanner, the reflectance value is preferred, which is the normalized and distance effect eliminated version of the intensity value.The coefficients <, >, H and I are constant for each scanner and the observations , J and M are variable for each point.
The equation part H = I is the linear distance error where H is the constant error and I is the proportional error whose contribution increases linearly as the distance increases.Coefficient d is a fractional number.
The function 9 : is the error contribution of the target reflectivity.It is a quadratic function of the distance .It is a piecewise function so that only the black (or absorbing) objects whose intensity I are less than the threshold : C contribute the error.
The linear distance error H = I and the target reflectivity error 9 : are the additive terms.The denominator term cos J is the error due to the incidence angle, which intensify or attenuate the additive error terms.
A single scan measurement configuration is proposed to derive the coefficients <, >, H and I. Two highly absorbent planar objects (black plates) and two perfectly reflective planar objects (white plates) were prepared.The pairs of black and white plates were located at close and far distances such as at 10 and 90 meters, respectively.Their orientations are kept perpendicular to the TLS so that the incidence angles become zero degrees.A single scan was performed.Assuming that the plates are exactly planar, the least squares plane fitting was computed for each plate using the conventional least squares parameter estimation method.The RMSE of the off-plane distances I was computed, represented with the symbol N.
The four N )T U , N )T C , N VT U and N VT C values are the RMS errors of the white plate at 10 meters, the black plate at 10 meters, the white plate at 90 meters and the black plate at 90 meters, respectively.Using the hypothetical graph in Figure 1, the scanner invariant parameters <, >, H and I can be computed straightforwardly.
Figure 1.Two sets of the black-and-white-plates are located at 10 and 90 meters distances, respectively.All of the four plates are at the rotation of 0 degrees of incidence angles.
Parameter H is the summation of the constant distance accuracy (W) provided by the manufacturer and N )T U which is the intercept of the red line (Figure 1).
The parameter I is the slope of the red line (Figure 1) which provides a kind of linear interpolation of the errors due to the distance.
The reflectivity error 9 : is computed only for the black objects.It is a kind of quadratic interpolation whose parameters < and > can be calculated by solving the following equation system: Once the scanner invariant coefficients <, >, H and I are computed, the range precision is calculated uniquely for each point with the Equations ( 5) and ( 6).Then, error ellipsoids are calculated (Figure 2).

THE LEAST ERRORS SURFACE RECONSTRUCTION
The proposed error ellipsoid is well defined mathematical representation of the physical reality, since it is formed by use of the most influential error sources.The size, shape and orientation of the error ellipsoid provide useful information about the metric quality of its point.
Every point is associated with its error ellipsoid whose semiaxis lengths are +( ) ,+( and+( * , respectively, where +( ) is the length of the semi-major axis.+( ) can be considered as a quality measure of the point.However, this approach is suboptimal, because +( ) is independent of the angular a priori precision values ( " and # ).The quality measure should be composed of both angular and the range precision values.The semi-diagonal axis length of the bounding box of the error ellipsoid (Figure 3) fulfils this requirement.It is computed as in the following equation: 11) In the TLS used projects, the target object has to be scanned from multiple viewpoints for full coverage.The consecutive scans should have the overlaps to some extent.In the overlapping parts, the same object surface is sampled at multiple times each of which are obviously at multiple quality levels.This situation results in data redundancy and the output surface mesh will not be in the uniform quality, if all the high and low quality points are used.
A decimation (subsampling) process should be performed in order that the points with the best spatial quality are selected over the redundant parts of the data.The remaining points are discarded.To achieve this objective, a voxel-based best quality point selection procedure was developed.The entire workflow is given in Figure 4.The workflow starts with the data acquisition and followed by the co-registration.In the data acquisition, stationary and repetitive 5 to maximum 10 scans are performed for the vertical and horizontal angle precision determination.In one of the scans, two sets of white and black plates, such as printed papers, are located at close and far distances for the range precision determination.These are the all extra labour costs for the point error model computation.
Figure 4.The processing workflow.
The scanner invariant vertical and horizontal angle precisions " and # , and the coefficients <, >, H and I are then computed.The every point in the every scan file is sought and its associated error ellipsoid is computed.All the scan files are merged, and a 3D voxelisation is established.The voxel sizes may vary depending on the object to be modelled and the required level of detail.The access and query mechanisms are provided by an octree data structure.The set of points lying in each voxel is determined.The number of points may vary for each voxel because the point density is not constant all over the point cloud.The final step is the selection of the point with the minimum a value for each voxel.As a result, the new point cloud is composed of points whose error ellipsoids have the minimum semi-diagonal axis length.The surface mesh is created using this decimated point cloud.

EXPERIMENTAL STUDY
The front façade of a building in the campus area of the Bulent Ecevit University (Zonguldak, Turkey) was scanned using a Faro Focus 3D X330 scanner.The building mostly contains planar components (Figure 5).Three scans were performed in such a way that all details of the façade can be covered.They were co-registered using Faro Scene v5.5, which is the bundled software of the Faro scanners.
The mean error of the co-registration was reported less than 2 mm by the software.The co-registration process was followed by computation of the error ellipsoids for each point of all three scans (Figure 6).Given the points and their error ellipsoids, the quality measure of each individual point was computed using Equation ( 11), which is the length of the semi-diagonal axis of the error ellipsoid.The computed point quality measures are visualized in Figure 7.
In Figure 7, the point quality is represented within a colour scale, varying between blue (min. 2 mm) and red (max.25 mm).
The blue colour corresponds to high quality and red to low quality points.Surface discontinuity, poor surface reflectance and higher incidence angles yield very low quality points.These points are represented by yellow-to-red colours.While in one scan the same part of the façade is measured with high quality points, in other scans the same part may be measured with low quality points.These points are mixed when the three scans are merged (Figure 8).The merged point cloud is messy and redundant, as it contains 5.2 million points.It is not optimal for surface mesh generation.A decimated point cloud was created using the voxel-based best point selection method.The voxel size was defined as 2 x 2 x 2 cm.The each voxel contributes only and only one point with the best spatial quality (Figure 9).When compared to the original merged point cloud (Figure 8), it is obvious that a great number of points of very low quality have been eliminated in the decimated point cloud (Figure 9), which has 1.4 million points.
However, there still exist some very low quality points shown in red colour which are particularly located at surface discontinuities.It is due to the fact that some voxels contain points all with very low quality, and even the best quality one is still in the low quality.
This situation requires an additional criterion.If the quality measure of the selected best quality point is greater than a threshold, i.e. if a d a , this point (and the associated voxel) is discarded.We empirically defined this threshold a as 6 mm.By threshold a value, we ensure that the surface mesh has a uniform positional quality of 6 mm or better at everywhere.
The newly decimated point cloud after the threshold a applied is shown in Figure 10.It has 1.1 million points now.
Figure 10.The final decimated point cloud.The points with quality measure of greater than 6 mm were discarded.
The final step of the workflow is the 3D surface mesh reconstruction.Geomagic Wrap 2015 and MeshLAB software were used to generate surface meshes both from the "original" (Figure 8) and "the best" (Figure 10) point clouds.Geomagic Wrap is a commercial software which can import point cloud data in various formats, and allows users to generate 3D surface meshes using a Delaunay-based algorithm (Edelsbrunner et al., 1998).On the other hand, MeshLab is a scientific software that can be used free of charge (Cignoni et al., 2008).Even though the MeshLab provides various mesh generation algorithms, we chose the ball-pivoting algorithm (Bernardini et al. 1999b) for this study.The results are shown in Figure 11 and Figure 12, respectively.Neither editing nor filtering operations have been applied on the results.
The "original" surface mesh (Figure 11a) has 9.8 million triangles, while "the best" surface mesh (Figure 11b) has only 1.9 million triangles.The proposed stochastic mesh generation method gives a better result with substantially less number of triangles by factor 5. The proposed method eliminates the noise significantly.Moreover, the redundant data is eliminated.Such a stochastic elimination of the erroneous points results in more accurate surface meshes.Identical results were obtained when using the academic software MeshLAB for the surface triangulation (Figure 12).zoom-in to the "original" mesh, (e-f) zoom-in to "the best" mesh.

CONCLUSIONS
We proposed a complete processing chain for high quality surface mesh generation specifically from the TLS derived point clouds.It is practical yet effective and powerful.The underlying methodology is generic, which can be implemented with using any scanner hardware and surface triangulation software combinations at no other operational cost.
The proposed point error model is anisotropic such that the spatial quality of each point is uniquely represented with its error ellipsoid.It is practical which can be computed in ordinary scanning projects with minimum labour and equipment costs (five more scans and four printed papers).The proposed surface reconstruction method is stochastic in which the spatial quality of the used points is rigorously regarded.The result is a nonredundant consensus surface with uniform positional quality.
Moreover, as the number of points is decreased, the computation time and memory requirements become feasible.
The proposed method can be favourably used in multi-sensor data fusion studies.
, & is the unit matrix, % = ( ) , ( , ( * contains the eigenvalues and ' = ) * is the eigenvector.The dimensions of the semi-axes of the error ellipsoid are the square The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2, 2018 ISPRS TC II Mid-term Symposium "Towards Photogrammetry 2020", 4-7 June 2018, Riva del Garda, Italy roots of the eigenvalues +( ) , +( , +( * .The axes orientation of the ellipsoid is given by the eigenvector.

Figure 2 .
Figure 2. Side view of a point cloud of two buildings which was acquired with a Riegl VZ400 scanner.The TLS stand point is labelled by the red cube.The estimated error ellipsoids are represented in green color.The ellipsoids are plotted at every 150-th points and their sizes are exaggerated for a better visualization.

Figure 3 .
Figure 3. Illustration of an error ellipsoid and the point quality metric a with respect to the TLS position.

Figure 5 .
Figure 5.The selected building façade for 3D surface mesh reconstruction.

Figure 6 .
Figure 6.Three point clouds of the building façade acquired by the Faro Focus 3D X330 scanner.The red cubes show the TLS stand points.The ellipsoids are plotted at every 200-th point and their sizes are exaggerated for a better visualization.

Figure 7 .
Figure 7. Visualization of the point quality of the three scans.The TLS stations are indicated by the white cubes.The same colour legend is used in Figure 8, 9 and 10.

Figure 8 .
Figure 8. Visualization of the point quality of the merged point cloud.

Figure 9 .
Figure 9. Visualization of the point quality of the decimated point cloud.The very low quality points (in red colour) still exist.

Figure 11 .Figure 12 .
Figure 11.3D surface meshes generated with the Geomagic Wrap 2015 software.(a) "Original" surface mesh, (b) "the best" surface mesh, (c-d) zoom-in to the "original" mesh, (e-f) zoom-in to "the best" mesh.