Mobile Laser Scanning Data for the Evaluation of Pavement Surface Distress

: The surface conditions of road pavements, including the occurrence and severity of distresses present on the surface, are an important indicator of pavement performance. Periodic monitoring and condition assessment is an essential requirement for the safety of vehicles moving on that road and the wellbeing of people. The traditional characterization of the di ﬀ erent types of distress often involves complex activities, sometimes ine ﬃ cient and risky, as they interfere with road tra ﬃ c. The mobile laser systems (MLS) are now widely used to acquire detailed information about the road surface in terms of a three-dimensional point cloud. Despite its increasing use, there are still no standards for the acquisition and processing of the data collected. The aim of our work was to develop a procedure for processing the data acquired by MLS, in order to identify the localized degradations that mostly a ﬀ ect safety. We have studied the data ﬂow and implemented several processing algorithms to identify and quantify a few types of distresses, namely potholes and swells / shoves, starting from very dense point clouds. We have implemented data processing in four steps: (i) editing of the point cloud to extract only the points belonging to the road surface, (ii) determination of the road roughness as deviation in height of every single point of the cloud with respect to the modeled road surface, (iii) segmentation of the distress (iv) computation of the main geometric parameters of the distress in order to classify it by severity levels. The results obtained by the proposed methodology are promising. The procedures implemented have made it possible to correctly segmented and identify the types of distress to be analyzed, in accordance with the on-site inspections. The tests carried out have shown that the choice of the values of some parameters to give as input to the software is not trivial: the choice of some of them is based on considerations related to the nature of the data, for others, it derives from the distress to be segmented. Due to the di ﬀ erent possible conﬁgurations of the various distresses it is better to choose these parameters according to the boundary conditions and not to impose default values. The test involved a 100-m long urban road segment, the surface of which was measured with an MLS installed on a vehicle that traveled the road at 10 km / h.


Introduction
The civil infrastructure pavement has a life cycle in many stages including design, construction (which includes new construction as well as preservation, maintenance, and rehabilitation), use and end-of-life. Since all pavements deteriorate over time, according to the Life Cycle Assessment (LCA) their monitoring is essential for road maintenance [1,2].
The reliable monitoring of the state of health of a road is an essential aspect of the infrastructure system, aimed at ensuring the safety and comfort of all users. The use of the most innovative systems jumps caused by road curbs. Balado et al. [52] introduced a method based on the use of a combination of geometrical and topological characteristics to classify some elements of street furniture (roads, sidewalks, treads, risers and curbs).
There are few extraction methods designed in situations where the edges of the road surface are not easily identifiable, mainly in rural and suburban areas. Among them, Yoon and Crane [53] identified the paved surface by computing the slopes and standard deviations of the elevation values of the single extracted profiles. Smadja et al. [35] use a second-degree polynomial regression to interpolate the points of a generic cross-section; Interpolation is made by applying the Random Sample Consensus (RANSAC) algorithm, which estimates the parameters of the mathematical model by iterative processes.
Recently, Yadav et al. [38] have proposed a method, designed for complex areas, divided into three steps: in the first step the elements and objects that are on the road surface and that do not belong to it are removed, thus dividing the ground points from the no-ground points; in the second step the topology and spatial density are analyzed, together with the intensity values to extract the paved surface; finally, in the last step the contours of the identified surface are analyzed and refined in order to obtain polylines of the pavement edges.
The aim of our work is to develop a procedure for processing the data acquired by MLS, starting from the filtering phase to extract only the points belonging to the ground up to the identification and quantification of some of the main localized degrades that most affect safety, namely swells, shoves and potholes.

Test Case and Data
The area to be surveyed is located in an urban environment in the municipality of Rome (Italy). The test was carried out on a stretch of road ( Figure 1). The road extends for about 100 m; it proceeds in a straight line for 80 m, then curves. The cross-section is typical of a two-lane road with a single carriageway 11 m wide; each lane is 3.50 m wide. Both sides of the road are used as parking.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 24 height jumps caused by road curbs. Balado et al. [52] introduced a method based on the use of a combination of geometrical and topological characteristics to classify some elements of street furniture (roads, sidewalks, treads, risers and curbs). There are few extraction methods designed in situations where the edges of the road surface are not easily identifiable, mainly in rural and suburban areas. Among them, Yoon and Crane [53] identified the paved surface by computing the slopes and standard deviations of the elevation values of the single extracted profiles. Smadja et al. [35] use a second-degree polynomial regression to interpolate the points of a generic cross-section; Interpolation is made by applying the Random Sample Consensus (RANSAC) algorithm, which estimates the parameters of the mathematical model by iterative processes.
Recently, Yadav et al. [38] have proposed a method, designed for complex areas, divided into three steps: in the first step the elements and objects that are on the road surface and that do not belong to it are removed, thus dividing the ground points from the no-ground points; in the second step the topology and spatial density are analyzed, together with the intensity values to extract the paved surface; finally, in the last step the contours of the identified surface are analyzed and refined in order to obtain polylines of the pavement edges.
The aim of our work is to develop a procedure for processing the data acquired by MLS, starting from the filtering phase to extract only the points belonging to the ground up to the identification and quantification of some of the main localized degrades that most affect safety, namely swells, shoves and potholes.

Test Case and Data
The area to be surveyed is located in an urban environment in the municipality of Rome (Italy). The test was carried out on a stretch of road ( Figure 1). The road extends for about 100 m; it proceeds in a straight line for 80 m, then curves. The cross-section is typical of a two-lane road with a single carriageway 11 m wide; each lane is 3.50 m wide. Both sides of the road are used as parking. The MLS data were acquired by a mobile laser scanner mounted on the roof of a vehicle driving on the road. The measuring head is equipped with two Riegl VQ-450 laser scanners, symmetrically configured on the left and right sides, pointing to the rear of the vehicle with an inclination angle of approximately 145°, as well as inertial and GNSS measurement equipment, housed under an aerodynamically-shaped protective cover. This configuration is called 'Butterfly' or 'X' pattern. Each The MLS data were acquired by a mobile laser scanner mounted on the roof of a vehicle driving on the road. The measuring head is equipped with two Riegl VQ-450 laser scanners, symmetrically configured on the left and right sides, pointing to the rear of the vehicle with an inclination angle of approximately 145 • , as well as inertial and GNSS measurement equipment, housed under an aerodynamically-shaped protective cover. This configuration is called 'Butterfly' or 'X' pattern. Each Remote Sens. 2020, 12, 942 5 of 25 VQ-450 scanner generates its own 360 • 'full circle' profile by the motorized scanning mechanism of the mirrors. The laser's acquisition frequency is 550 kHz. The main characteristics of the MLS used are listed in Table 1. A complete road survey was carried out at a driving speed of about 10 km/h. The output data are point clouds directly georeferenced in the UTM/ETRS00 cartographic system [54].
The density of the points of the MLS cloud depends both on the instrument-target center distance and on the angle of incidence. For an average distance of 5 m, at a vehicle speed of 10 km/h, the value of the point-density is estimated at 32,544 points/m 2 , while it decreases to 4688 for a speed of 60 km/h. Figure 2 shows the actual density maps (not estimated). The distance between the scanning lines is approximately 1.2 cm and 7.0 cm, respectively for driving speeds of 10 and 60 km/h. Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 24 VQ-450 scanner generates its own 360° 'full circle' profile by the motorized scanning mechanism of the mirrors. The laser's acquisition frequency is 550 kHz. The main characteristics of the MLS used are listed in Table 1.  [54].
The density of the points of the MLS cloud depends both on the instrument-target center distance and on the angle of incidence. For an average distance of 5 m, at a vehicle speed of 10 km/h, the value of the point-density is estimated at 32,544 points/m 2 , while it decreases to 4688 for a speed of 60 km/h. Figure 2 shows the actual density maps (not estimated). The distance between the scanning lines is approximately 1.2 cm and 7.0 cm, respectively for driving speeds of 10 and 60 km/h. Follow-up analyses were made on the densest cloud, the one acquired in only one travel direction, at a driving speed of 10 km/h.

Methods
The proposed method aims both at identifying the relevant distresses and at estimating their size.
Data processing has been carried out in four main steps ( Figure 3 shows a workflow of the procedure):  the first step focuses on editing the MLS point cloud to extract only the points belonging to the road surface;  the second step consists in evaluating the roughness, determining the height deviation of every single point of the cloud with respect to the reconstructed model of the pavement surface; Follow-up analyses were made on the densest cloud, the one acquired in only one travel direction, at a driving speed of 10 km/h.

Methods
The proposed method aims both at identifying the relevant distresses and at estimating their size. Data processing has been carried out in four main steps ( Figure 3 shows a workflow of the procedure): • the first step focuses on editing the MLS point cloud to extract only the points belonging to the road surface; • the second step consists in evaluating the roughness, determining the height deviation of every single point of the cloud with respect to the reconstructed model of the pavement surface; • the third step involves the creation of a binary image starting from the built DEM; two different images will be created and analyzed separately, one for positive displacements, the other for negative ones; Remote Sens. 2020, 12, 942 6 of 25 • the last step, the fourth, focuses on the computation of the main geometric parameters of each segmented region. According to these, each region will be classified according to the severity levels given in the standards.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 24  the third step involves the creation of a binary image starting from the built DEM; two different images will be created and analyzed separately, one for positive displacements, the other for negative ones;  the last step, the fourth, focuses on the computation of the main geometric parameters of each segmented region. According to these, each region will be classified according to the severity levels given in the standards.

Semi-automatic Extraction of the Road Surface
The first step of the whole process that will lead to the identification and classification of the distress is the editing of the point cloud to remove the outliers and extract the bare road surface, as shown in the workflow of Figure 3. This phase consists of four steps: 1. Extracting the vehicle's trajectory from the LAS data To extract from the whole cloud of laser points only a few points that will be the vertexes of the polyline representative of the vehicle's trajectory, the points acquired with a given scanning angle are extracted. The matrix containing the planimetric coordinates of the points acquired with a scanning angle equal to 0° has been ordered, in ascending order, according to GPS time (value measured for each single point). The first vertex V1 of the polyline is identified by the planimetric coordinates contained in the first line of the matrix thus ordered. Starting from the first vertex V1 all points contained in a sphere with center at the point itself and radius equal to the length of the polyline segments are extracted. The function used is the 'findNeighborsInRadius'. In our test case, the polyline segments are 5 m long. The function outputs the line indices of the matrix containing the planimetric coordinates of the points inside the sphere and the distances of each point from its center. The next vertex V2 of the polyline is identified by the planimetric coordinates of the point having the greatest distance from V1 and the largest GPS time. V2 is the center of the next research sphere. The cycle repeats until the coordinates of all polyline vertices are found. Figure 4 shows a schematic drawing of the procedure for constructing the polyline that discretizes the vehicle's trajectory.

Semi-Automatic Extraction of the Road Surface
The first step of the whole process that will lead to the identification and classification of the distress is the editing of the point cloud to remove the outliers and extract the bare road surface, as shown in the workflow of Figure 3. This phase consists of four steps:

Extracting the vehicle's trajectory from the LAS data
To extract from the whole cloud of laser points only a few points that will be the vertexes of the polyline representative of the vehicle's trajectory, the points acquired with a given scanning angle are extracted. The matrix containing the planimetric coordinates of the points acquired with a scanning angle equal to 0 • has been ordered, in ascending order, according to GPS time (value measured for each single point). The first vertex V 1 of the polyline is identified by the planimetric coordinates contained in the first line of the matrix thus ordered. Starting from the first vertex V 1 all points contained in a sphere with center at the point itself and radius equal to the length of the polyline segments are extracted. The function used is the 'findNeighborsInRadius'. In our test case, the polyline segments are 5 m long. The function outputs the line indices of the matrix containing the planimetric coordinates of the points inside the sphere and the distances of each point from its center. The next vertex V 2 of the polyline is identified by the planimetric coordinates of the point having the greatest distance from V 1 and the largest GPS time. V 2 is the center of the next research sphere. The cycle repeats until the coordinates of all polyline vertices are found. Figure 4 shows a schematic drawing of the procedure for constructing the polyline that discretizes the vehicle's trajectory.

Extracting sub-clouds inside built-up boxes
To reduce the CPU time, the whole point cloud has been divided into many sub-clouds, whose points are contained in polygonal boxes purposely built to divide it. The vertices of the boxes are aligned in a direction orthogonal to the polyline that describes the trajectory of the vehicle.
Given the n vertices V of the polyline, whose coordinates are contained in the matrix where im is the imaginary unit, i.e., im = 0.0000 + 1.0000•i.
Then, the row vector , which contains the coordinate differences between the vth vertex and the vth-1, is built: Then it is built the row vector where M1 indicates the first element of M, M1:n-2 indicates the elements of M from the 1 st to the n th -2 (with n = number of vertices) and Mn-1 indicates the last element of the row vector M, which has the dimensions (n-1). For each pair of planimetric coordinates Ev, Nv of the generic vertex v the two pairs of offset where:

Extracting sub-clouds inside built-up boxes
To reduce the CPU time, the whole point cloud has been divided into many sub-clouds, whose points are contained in polygonal boxes purposely built to divide it. The vertices of the boxes are aligned in a direction orthogonal to the polyline that describes the trajectory of the vehicle.
Given the n vertices V of the polyline, whose coordinates are contained in the matrix V ∈ R n,2 , to compute the coordinates of the offset points, a row vector M ∈ R 1,n is built, which contains the coordinates E v , N v (East and North respectively) of each vertex v of the polyline: where im is the imaginary unit, i.e., im = 0.0000 + 1.0000·i. Then, the row vector ∆M ∈ R 1,n−1 , which contains the coordinate differences between the v th vertex and the v th−1 , is built: Then it is built the row vector ∆M II ∈ R 1,n : where ∆M 1 indicates the first element of ∆M, ∆M 1:n−2 indicates the elements of ∆M from the 1st to the nth−2 (with n = number of vertices) and ∆M n−1 indicates the last element of the row vector ∆M, which has the dimensions (n − 1). For each pair of planimetric coordinates E v , N v of the generic vertex v the two pairs of offset coordinates E o ∈ R 1,2 , N o ∈ R 1,2 are obtained by subtracting and adding the scalar of the real or imaginary part of the respective row of the matrix O ∈ R 1,n .
Once a certain offset value d is set, the coordinates of the vertices of the box are: where: Remote Sens. 2020, 12, 942 8 of 25 The 'offset points' correspond to the vertices of the quadrilateral (the boxes). The length of the boxes is that of the segments of the polyline, the width varies depending on the width of the carriageway. In our case, the dimensions of the boxes are 10 × 5 m. Figure 5 shows a scheme of the construction of the boxes.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 24 The 'offset points' correspond to the vertices of the quadrilateral (the boxes). The length of the boxes is that of the segments of the polyline, the width varies depending on the width of the carriageway. In our case, the dimensions of the boxes are 10 × 5 m. Figure 5 shows a scheme of the construction of the boxes.

Down-sampling of the point cloud
The next step consists of down-sampling the point cloud to obtain the nodes that are used as the center of the research spheres for the follow-up road surface extraction process. To do this, the Matlab function 'pcdownsample', method 'gridAverage', has been used. A regular voxel grid is built with a user assigned step. The height values of the points surrounding the grid node are averaged and the average is assigned as the height value to the node. Before proceeding with the editing process, nodes in areas outside the road surface must be removed. By using the 'pcnormals' function, for each node in the down-sampled grid, a plane has been fitted by least squares and the normal vector to it has been estimated. Points whose normal vector has a component along the z axis (with z in the direction of the ellipsoidal height) of less than 1 (with a margin of error of 5 × 10 −3 ) have been removed. Finally, the 'pcdenoise' function has been applied to remove isolated or edge nodes.

Road surface extraction
In the last editing step, we used the M-estimator SAmple Consensus (MSAC) algorithm, which is a robust variant of the Random Sample Consensus (RANSAC) algorithm, an iterative method to estimate the parameters of a mathematical model from a set of observed data that contains outliers [55].
The processing steps are:  Starting from the remaining nodes of the voxel grid, using the 'findNeighborsInRadius' function, a sphere centered in each node is built. The radius of the sphere has been chosen to be proportional to the grid step; it is recommended to choose a radius value such that two adjacent spheres have a large overlapping area, for example, three times the grid step.  Using the 'pcfitplane' function, the points (of the original point cloud) falling within the sphere are fitted on a plane, the parameters of whose equation are estimated. The variables of the function are: the 'maxDistance' between the plane and the generic point so that it can be considered an inlier, the 'referenceVector' which is a constraint on the orientation of the

Down-sampling of the point cloud
The next step consists of down-sampling the point cloud to obtain the nodes that are used as the center of the research spheres for the follow-up road surface extraction process. To do this, the Matlab function 'pcdownsample', method 'gridAverage', has been used. A regular voxel grid is built with a user assigned step. The height values of the points surrounding the grid node are averaged and the average is assigned as the height value to the node. Before proceeding with the editing process, nodes in areas outside the road surface must be removed. By using the 'pcnormals' function, for each node in the down-sampled grid, a plane has been fitted by least squares and the normal vector to it has been estimated. Points whose normal vector has a component along the z axis (with z in the direction of the ellipsoidal height) of less than 1 (with a margin of error of 5 × 10 −3 ) have been removed. Finally, the 'pcdenoise' function has been applied to remove isolated or edge nodes.

Road surface extraction
In the last editing step, we used the M-estimator SAmple Consensus (MSAC) algorithm, which is a robust variant of the Random Sample Consensus (RANSAC) algorithm, an iterative method to estimate the parameters of a mathematical model from a set of observed data that contains outliers [55].
The processing steps are: • Starting from the remaining nodes of the voxel grid, using the 'findNeighborsInRadius' function, a sphere centered in each node is built. The radius of the sphere has been chosen to be proportional to the grid step; it is recommended to choose a radius value such that two adjacent spheres have a large overlapping area, for example, three times the grid step. • Using the 'pcfitplane' function, the points (of the original point cloud) falling within the sphere are fitted on a plane, the parameters of whose equation are estimated. The variables of the function are: the 'maxDistance' between the plane and the generic point so that it can be considered an inlier, the 'referenceVector' which is a constraint on the orientation of the reference plane, and the 'maxAngularDistance', which is a threshold value for the angular distance between the normal vector of the fitted plane and another vector (in the direction of the ellipsoidal height). Since the plane has been created with the iterative MSAC process, its position cannot be fixed a priori, the plane must be able to rotate in space: this is why it is required to set a threshold value for the 'maxAngularDistance'. Moreover, since in our case we want to identify surfaces with mostly "horizontal" position, a 'referenceVector' dim [0,0,1] has been used.

•
Once the outliers have been removed, the cycle starts again on one other grid node until the editing process is complete.

Roughness Evaluation
The evaluation of localized distress in 3D space is rather complex; there is no ideal reference surface to model so that any anomalies present can be measured with respect to it. It could be considered to use a plane to characterize a single lane and determine the height variations of the road surface with respect to this one [43]. Nevertheless, there are essentially two main problems in this approach: it is not trivial to choose the longitudinal dimensions that the reference plane should assume; the sides of a plane must coincide with the sides of adjacent planes; this is a very rigid constraint.
That is why we have decided to use several local reference planes, built gradually for each point of the cloud belonging to the road surface. The procedure follows the following steps: • A sphere with an assigned radius (kernel size) and center at each point of the cloud has been built ( Figure 6).

•
All the points inside the sphere are fitted on a plane. The 'findNeighborsInRadius' function has been used to select the points. The plane has been fitted to data with a least-squares robust regression, the 'Bisquare Weights'. This method minimizes a weighted sum of squares, where the weight of each data point depends on its distance from the fitted plane; the farther away is the point, the less weight it gets. • Finally, the roughness value ∆h is computed as the distance between the point (center of the sphere) and the fitting plane. Positive values of ∆h indicate that the point is below the plane (downward deformations) while negative values are defining the upward deformations.
Besides the computation of ∆h, the algorithm provides in output, as an estimation of the goodness of the fitting, the RMSE (Root Mean Square Error) value, that is also associated to all the point of the cloud; its layout contributes to give a judgment of the road surface conditions. Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 24 reference plane, and the 'maxAngularDistance', which is a threshold value for the angular distance between the normal vector of the fitted plane and another vector (in the direction of the ellipsoidal height). Since the plane has been created with the iterative MSAC process, its position cannot be fixed a priori, the plane must be able to rotate in space: this is why it is required to set a threshold value for the 'maxAngularDistance'. Moreover, since in our case we want to identify surfaces with mostly "horizontal" position, a 'referenceVector' dim [0,0,1] has been used.  Once the outliers have been removed, the cycle starts again on one other grid node until the editing process is complete.

Roughness Evaluation
The evaluation of localized distress in 3D space is rather complex; there is no ideal reference surface to model so that any anomalies present can be measured with respect to it. It could be considered to use a plane to characterize a single lane and determine the height variations of the road surface with respect to this one [43]. Nevertheless, there are essentially two main problems in this approach: it is not trivial to choose the longitudinal dimensions that the reference plane should assume; the sides of a plane must coincide with the sides of adjacent planes; this is a very rigid constraint.
That is why we have decided to use several local reference planes, built gradually for each point of the cloud belonging to the road surface. The procedure follows the following steps:


A sphere with an assigned radius (kernel size) and center at each point of the cloud has been built ( Figure 6).  All the points inside the sphere are fitted on a plane. The 'findNeighborsInRadius' function has been used to select the points. The plane has been fitted to data with a least-squares robust regression, the 'Bisquare Weights'. This method minimizes a weighted sum of squares, where the weight of each data point depends on its distance from the fitted plane; the farther away is the point, the less weight it gets.  Finally, the roughness value h is computed as the distance between the point (center of the sphere) and the fitting plane. Positive values of h indicate that the point is below the plane (downward deformations) while negative values are defining the upward deformations. Besides the computation of h, the algorithm provides in output, as an estimation of the goodness of the fitting, the RMSE (Root Mean Square Error) value, that is also associated to all the point of the cloud; its layout contributes to give a judgment of the road surface conditions.

Distress Identification and Evaluation
The quantification and classification of the different distress of the road surface has been done using image segmentation algorithms, treating the point cloud as a digital image in which pixels are grouped according to specific conditions. The implemented procedure follows the following steps: 1. Interpolation of the roughness values h on the nodes of a grid

Distress Identification and Evaluation
The quantification and classification of the different distress of the road surface has been done using image segmentation algorithms, treating the point cloud as a digital image in which pixels are grouped according to specific conditions. The implemented procedure follows the following steps: 1. Interpolation of the roughness values ∆h on the nodes of a grid The coordinates of the point cloud (N, E, ∆h) have been interpolated to build a raster DEM where the coordinate h is replaced with the roughness value ∆h. For the building of the DEM the IDW (Inverse Distance Weighted) interpolator (with distance power 2) implemented in Surfer ® ver. 12 by Golden Software was used [56]. We chose it because, although it is a very common deterministic interpolator, it works well for high-density LiDAR data [57]; moreover, in [58] some comparisons made with reference DEM have shown a better correlation when using IDW than using other more robust interpolation methods.
The spatial resolution ρ of the DEM was chosen according to two parameters that characterize any scan: the minimum surface density value and the mean shortest distance between the scan lines. The mathematical formula that gives the resolution as a function of surface density, for a quite regular distribution of points such as the test-case is [59]: where A is the surface area of the scanned site and N is the total number of points in the area A. The simplified formula that gives the grid density as a function of the average distance between adjacent scan lines is [59]: where d i,j is the average distance between scan lines. The resolution of the DEM must be equal to or greater than the maximum value between those computed in this way. In our test case the (6) gives 0.8 cm in output and the (7) one gives 1.1 cm. The resolution of the DEM we have chosen is 1.5 cm.

Creating a binary image
The raster DEM, in which the pixel size corresponds to the grid step, has been transformed into a binary image. To the pixels with ∆h greater in absolute value of a given threshold value (different for upward and downward displacements) is assigned the value 1 (that corresponds to white) while to all the other pixels is assigned the value 0 (black pixels). The white pixels will, therefore, be those affected by distress while the black pixels will be those not affected. Two separate binary images are created for positive and negative roughness values. Figure 7a shows an example of a section of road affected by downward displacements.
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 24 The coordinates of the point cloud (N, E, Δh) have been interpolated to build a raster DEM where the coordinate h is replaced with the roughness value Δh. For the building of the DEM the IDW (Inverse Distance Weighted) interpolator (with distance power 2) implemented in Surfer® ver. 12 by Golden Software was used [56]. We chose it because, although it is a very common deterministic interpolator, it works well for high-density LiDAR data [57]; moreover, in [58] some comparisons made with reference DEM have shown a better correlation when using IDW than using other more robust interpolation methods.
The spatial resolution  of the DEM was chosen according to two parameters that characterize any scan: the minimum surface density value and the mean shortest distance between the scan lines. The mathematical formula that gives the resolution as a function of surface density, for a quite regular distribution of points such as the test-case is [59]: where A is the surface area of the scanned site and N is the total number of points in the area A. The simplified formula that gives the grid density as a function of the average distance between adjacent scan lines is [59]: where di,j is the average distance between scan lines. The resolution of the DEM must be equal to or greater than the maximum value between those computed in this way. In our test case the (6) gives 0.8 cm in output and the (7) one gives 1.1 cm. The resolution of the DEM we have chosen is 1.5 cm.

Creating a binary image
The raster DEM, in which the pixel size corresponds to the grid step, has been transformed into a binary image. To the pixels with h greater in absolute value of a given threshold value (different for upward and downward displacements) is assigned the value 1 (that corresponds to white) while to all the other pixels is assigned the value 0 (black pixels). The white pixels will, therefore, be those affected by distress while the black pixels will be those not affected. Two separate binary images are created for positive and negative roughness values. Figure 7a shows an example of a section of road affected by downward displacements.

Computation of geometric parameters and segmentation
The next step is the implementation of both geometric and morphological operators in order to remove or group isolated white pixel regions with certain characteristics; this is the proper segmentation. The main geometric parameters computed for each region of the binary matrix are: perimeter, area, volume (Figure 8).

Computation of geometric parameters and segmentation
The next step is the implementation of both geometric and morphological operators in order to remove or group isolated white pixel regions with certain characteristics; this is the proper segmentation. The main geometric parameters computed for each region of the binary matrix are: perimeter, area, volume ( Figure 8). Given the coordinates (xp,k,yp,k) of the vertices of the n segments of the polyline that makes the outline of the distress, the perimeter is computed as: The area has been computed as the sum of the pixels of each region. The volume has been computed as: where ij is the pixel size, hj is the height of each pixel and npx is the total number of pixels in each region.
In this step, other geometrical parameters were also computed, such as the length of the axes of the circumscribed ellipse ('regionprops' function) to be used together with the previous ones for the identification of the zones to be segmented. Their selection is based on the values assigned to the parameters set in accordance with the specifications provided by ASTM-D6433 [11].
In particular, the following parameters have been analyzed for the segmentation:  h values. The regions will be divided according to the level of severity, which in turn depends on the depths or heights of the distress.  Area. The level of severity of distress is a function not only of its dimension in elevation but also of its extension in planimetry.


Length of the axes of the circumscribed ellipse and the mean diameter ('regionprops' function). In addition to the extent that distress may have, the other parameters that affect the severity levels are the mean diameter and orientation with respect to the direction of travel. The geometrical parameters mentioned above have been combined and analyzed according to the type of distress (downward or upward displacements) that is to be characterized. The removal of regions that do not meet the requirements is based on a cycle of conditional instructions, implemented using the 'RemoveMask' function.
Finally, expansion and erosion operators have been applied on the binary image to add or remove pixels on the boundaries of the regions in order to obtain a regular contour of the distress. Contours of the regions have been identified using the 'bwboundaries' function, based on the 'Moore-Neighbor' tracking algorithm, modified by Jacob's stopping criterion [60] and the morphological functions 'strel' and 'imfill' have been implemented [61,62]. Figure 7b shows an example of the output.

Extraction of the Road Surface
This phase aims at removing the outliers and extracting only the surface of the road pavement. For computational constraints, the whole point cloud has been divided into 25 polygonal boxes each with a size of 12 × 5 m, built around the trajectory of the vehicle reconstructed by extracting from the Given the coordinates (x p,k ,y p,k ) of the vertices of the n segments of the polyline that makes the outline of the distress, the perimeter is computed as: The area has been computed as the sum of the pixels of each region. The volume has been computed as: where ∆i j is the pixel size, ∆h j is the height of each pixel and npx is the total number of pixels in each region. In this step, other geometrical parameters were also computed, such as the length of the axes of the circumscribed ellipse ('regionprops' function) to be used together with the previous ones for the identification of the zones to be segmented. Their selection is based on the values assigned to the parameters set in accordance with the specifications provided by ASTM-D6433 [11].
In particular, the following parameters have been analyzed for the segmentation: • ∆h values. The regions will be divided according to the level of severity, which in turn depends on the depths or heights of the distress.

•
Area. The level of severity of distress is a function not only of its dimension in elevation but also of its extension in planimetry.

•
Length of the axes of the circumscribed ellipse and the mean diameter ('regionprops' function). In addition to the extent that distress may have, the other parameters that affect the severity levels are the mean diameter and orientation with respect to the direction of travel.
The geometrical parameters mentioned above have been combined and analyzed according to the type of distress (downward or upward displacements) that is to be characterized. The removal of regions that do not meet the requirements is based on a cycle of conditional instructions, implemented using the 'RemoveMask' function.
Finally, expansion and erosion operators have been applied on the binary image to add or remove pixels on the boundaries of the regions in order to obtain a regular contour of the distress. Contours of the regions have been identified using the 'bwboundaries' function, based on the 'Moore-Neighbor' tracking algorithm, modified by Jacob's stopping criterion [60] and the morphological functions 'strel' and 'imfill' have been implemented [61,62]. Figure 7b shows an example of the output.

Extraction of the Road Surface
This phase aims at removing the outliers and extracting only the surface of the road pavement. For computational constraints, the whole point cloud has been divided into 25 polygonal boxes each with a size of 12 × 5 m, built around the trajectory of the vehicle reconstructed by extracting from the MLS point cloud the only points with a scanning angle of 0 degree (with a margin of error of 5 × 10 −3 ). These points were interpolated linearly to construct a polyline made of 5 m long sections. Only 20 boxes, (from 5th to 24th) were considered for the next editing step. Figure 9 shows the 25 polygonal boxes built to divide the cloud into sub-clouds (a); in (b) a zoomed excerpt of the points of the cloud corresponding to the MLS track is shown.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 24 MLS point cloud the only points with a scanning angle of 0 degree (with a margin of error of 5 × 10 −3 ). These points were interpolated linearly to construct a polyline made of 5 m long sections. Only 20 boxes, (from 5th to 24th) were considered for the next editing step. Figure 9 shows the 25 polygonal boxes built to divide the cloud into sub-clouds (a); in (b) a zoomed excerpt of the points of the cloud corresponding to the MLS track is shown. On the points belonging to each box, a voxel grid (step size 0.20 m) has been built in order to down-sample the original point cloud. On the nodes of the voxel grid the RANSAC-based plane fitting algorithm was run, and a sphere with a diameter 3 times the grid step, namely 0.60 m, has been built to remove the outliers.
The following were considered as outliers:  the points that have a distance from the plane greater than the selected threshold. In our test we set a threshold absolute value of 2.5 cm, within the range of values corresponding to the pavement's macro-texture;  points belonging to a plane whose normal has a z component that differs angularly more than the selected threshold from the direction of the ellipsoidal height. To choose the input value for this parameter, namely 'maxAngularDistance', we made a few tests. Given unchanged values for the other parameters, most damaged areas were removed for angles less than 5°. For angles greater than 10°, some points not belonging to the road surface were not removed, in detail along the edges and in correspondence of objects present on the road surface. In our test case, the value that led to better results was 10°. Figure 10 shows a box, number 3, in profile and in axonometric view, before and after editing. Figure 11 shows the whole point cloud, either original or edited.  On the points belonging to each box, a voxel grid (step size 0.20 m) has been built in order to down-sample the original point cloud. On the nodes of the voxel grid the RANSAC-based plane fitting algorithm was run, and a sphere with a diameter 3 times the grid step, namely 0.60 m, has been built to remove the outliers.
The following were considered as outliers: • the points that have a distance from the plane greater than the selected threshold. In our test we set a threshold absolute value of 2.5 cm, within the range of values corresponding to the pavement's macro-texture; • points belonging to a plane whose normal has a z component that differs angularly more than the selected threshold from the direction of the ellipsoidal height. To choose the input value for this parameter, namely 'maxAngularDistance', we made a few tests. Given unchanged values for the other parameters, most damaged areas were removed for angles less than 5 • . For angles greater than 10 • , some points not belonging to the road surface were not removed, in detail along the edges and in correspondence of objects present on the road surface. In our test case, the value that led to better results was 10 • . Figure 10 shows a box, number 3, in profile and in axonometric view, before and after editing. Figure 11 shows the whole point cloud, either original or edited. value for this parameter, namely 'maxAngularDistance', we made a few tests. Given unchanged values for the other parameters, most damaged areas were removed for angles less than 5°. For angles greater than 10°, some points not belonging to the road surface were not removed, in detail along the edges and in correspondence of objects present on the road surface. In our test case, the value that led to better results was 10°. Figure 10 shows a box, number 3, in profile and in axonometric view, before and after editing. Figure 11 shows the whole point cloud, either original or edited.

Roughness Evaluation
The test was carried out on an area characterized by a variety of degradation, at different levels of severity, chosen to test and validate the proposed method. In this phase was computed, as described in Section 3.2, the difference in height of each single point belonging to the point cloud with respect to the reference plane, constructed by interpolating on it the points belonging to a sphere of a given radius. The only parameter to set is, therefore, the radius of the sphere (the size of the kernel). The kernel size has been set to 0.60 m, congruent with some distress evaluation.
The result is a point cloud where each point is associated with a scalar representing the distance between the point and the corresponding plane. Positive values indicate downward displacements (potholes) while negative values indicate upward displacements (swells/shoves). Figure 12 shows the roughness value h computed for each point of the cloud, the grey tone represents its value.
Positive values (towards white) represent the points below the reference plane, while negative values (towards black) represent the points above it. The maximum value found for the potholes was 46 mm, while for the swells/shoves was 41 mm.

Roughness Evaluation
The test was carried out on an area characterized by a variety of degradation, at different levels of severity, chosen to test and validate the proposed method. In this phase was computed, as described in Section 3.2, the difference in height of each single point belonging to the point cloud with respect to the reference plane, constructed by interpolating on it the points belonging to a sphere of a given radius. The only parameter to set is, therefore, the radius of the sphere (the size of the kernel). The kernel size has been set to 0.60 m, congruent with some distress evaluation.
The result is a point cloud where each point is associated with a scalar representing the distance between the point and the corresponding plane. Positive values indicate downward displacements (potholes) while negative values indicate upward displacements (swells/shoves). Figure 12 shows the roughness value ∆h computed for each point of the cloud, the grey tone represents its value.
Positive values (towards white) represent the points below the reference plane, while negative values (towards black) represent the points above it. The maximum value found for the potholes was 46 mm, while for the swells/shoves was 41 mm.

Distress Identification
The point cloud with the roughness values put in place of the heights was interpolated on a DEM grid using the software package Surfer ® ver. 12 by Golden software [56]. The inverse distance interpolation method (power 2) of the grid was used. The grid step was set to 1.5 cm, depending on the density of the MLS data. The search radius was set to 0.50 m. The output grid is a matrix, in which a roughness value is assigned to each cell. The next step concerns the segmentation for the extraction and cataloging of surface distress, according to the severity levels introduced by the standard [11]. The estimated distresses have been divided into two main categories:  downward displacements, that includes depressions, cracking, potholes, rutting, etc.  upward displacements, that includes slides, bumps, swells, shoves, patches, etc.   Figure 13 shows the Root Mean Square Error (RMSE) associated with the points in the cloud. It was computed to assess the goodness of fit, it is also significant for an initial classification of the pavement condition. Areas with low RMSE values should correspond to regular paving conditions (good road surface condition), while high values show surface irregularities (poor condition).

Distress Identification
The point cloud with the roughness values put in place of the heights was interpolated on a DEM grid using the software package Surfer ® ver. 12 by Golden software [56]. The inverse distance interpolation method (power 2) of the grid was used. The grid step was set to 1.5 cm, depending on the density of the MLS data. The search radius was set to 0.50 m. The output grid is a matrix, in which a roughness value is assigned to each cell. The next step concerns the segmentation for the extraction and cataloging of surface distress, according to the severity levels introduced by the standard [11]. The estimated distresses have been divided into two main categories:  downward displacements, that includes depressions, cracking, potholes, rutting, etc.  upward displacements, that includes slides, bumps, swells, shoves, patches, etc.

Distress Identification
The point cloud with the roughness values put in place of the heights was interpolated on a DEM grid using the software package Surfer ® ver. 12 by Golden software [56]. The inverse distance interpolation method (power 2) of the grid was used. The grid step was set to 1.5 cm, depending on the density of the MLS data. The search radius was set to 0.50 m. The output grid is a matrix, in which a roughness value is assigned to each cell. The next step concerns the segmentation for the extraction and cataloging of surface distress, according to the severity levels introduced by the standard [11]. The estimated distresses have been divided into two main categories: • downward displacements, that includes depressions, cracking, potholes, rutting, etc. • upward displacements, that includes slides, bumps, swells, shoves, patches, etc.
Among the downward displacements we have focused our analysis on the potholes, which are small, usually less than 750 mm in diameter, bowl-shaped holes in the pavement surface. To decide the parameters to be set for segmentation, we have assumed that they are not to be defined as such if they have a diameter less than 100 mm and severity levels are significant when the depth exceeds 13 mm. In addition, standards and literature define them as dangerous when they are more than 1 inch (about 25 mm) deep and the area is greater than 1 square foot (about 1 dm 2 , which corresponds to an average diameter of about 55 cm) [63].
Consequently, only regions with axes smaller than 100 mm, area smaller than 1 dm 2 and roughness values ∆h, positive, smaller than 13 mm have been discarded.
Among the upward displacements, we have focused our analysis on swells and shoves. Swell is characterized by un upward bulge in the pavement's surface whereas shove is a permanent, longitudinal displacement/ripple of a localized area of the pavement surface. Swelling looks like as a long, gradual waves more than 3 m long and shows as short, abrupt waves on the pavement surface. Both can make the ride uncomfortable and unsafe, producing significant damage in the long term. These types of distress tend to be mostly functional, i.e., the superstructure is still efficient but has some critical issues in terms of regularity and adherence.
As far as swells/shoves are concerned, only regions with an area smaller than 10 dm 2 and |∆h| larger than 5 mm (∆h negative) were discarded. Figure 14 shows the identified distresses, which are 29 potholes and 98 swells/shoves.
In the test case, the clearly recognizable swells are the number 25, 35, 38, 45, 49, 66 and 72 (Figure 14b) mainly oriented along the direction of travel and longer than 3 m. The shoves are all isolated segmented regions. Most were formed by the longitudinal acceleration (braking, speed up) of the tire and the vehicle load. Some of these, on the other hand, was formed due to poor compaction and the thrust (by the vehicles) of the superficial part of the pavement against the hydraulic devices consisting of more rigid elements (16, 57, 77, in Figure 14b), their shapes clearly outline manholes. The segmented regions are congruent with the particular configuration of the type of distress, i.e., swollen and isolated areas with longitudinal development, as shown in Figure 14b, numbers 10, 24, 22, 19, 60 and 73.
In addition, the distresses have been classified according to the severity levels deduced from ASTM D6433 [11], Table 2. The severity levels are divided into low (L), medium (M) and high (H). By adopting a speed of travel suited to the type of road being examined, a low level of severity corresponds to a situation where a reduction in speed is not necessary to respect good levels of comfort and safety, even with the presence of different types of distress. Localized distresses can create vibrations and jerks of the vehicle that only cause minor discomforts to the driver. A medium level is synonymous with significant vehicle vibration; speed reduction is necessary to ensure good levels of safety and comfort. A high level of severity highlights a situation where vibrations and jerks are so strong that the user is forced to significantly reduce the transit speed in order to achieve an acceptable level of comfort and safety; localized distresses create significant discomfort, as well as a high risk of damage to the vehicle.
In the test case, the most concerned are 21 potholes, in particular, 16 with a medium severity level and five with a high severity level. As far as swells/shoves are concerned, most (93) were found to be of low severity level and only five of medium severity level. During maintenance, in addition to the category of road that also determines the degree of priority, the geometric characteristics of the distress provide the parameters used to estimate the intervention time limit. The maximum time allowed between identification and repair depends on two factors: the size of the pothole and the hierarchy of the road. Both factors have an impact on the hazard and the level of risk created by the pothole itself. Potholes with a high level of severity, for main roads (e.g., motorways) should be repaired within 5 h after identification [63]. At the same time, potholes with a 'medium' level acquire different weights as the level of the hierarchy of the road changes.
For all identified distresses the main geometric parameters have been computed, i.e., for the potholes the perimeter, area, volume, maximum depth, axes of the circumscribed ellipse and the mean diameter, while for the swells/shoves only the perimeter, area, volume and maximum depth. Table 3 shows, for some significant examples of identified distresses (both potholes and swells/shoves), the computed geometric parameters.
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 24 transit speed in order to achieve an acceptable level of comfort and safety; localized distresses create significant discomfort, as well as a high risk of damage to the vehicle.

Discussion
The results obtained with the proposed method depend mainly on the value of the parameters given as input to the software. In particular, in the editing phase, it is important to choose the correct values because the reconstructed 3D model of the road surface represents the reference surface on which the distresses will be highlighted.
The choice of some parameters is based on considerations related to the nature of the input data: the DEM resolution, for example, is in relation to the characteristics of the point cloud, i.e., its density and accuracy. Other parameters are decided according to the type of distress to be identified. For others, some tests have been carried out to verify the influence of the choice of certain values on the results. In the segmentation procedure (see Section 3.3), the parameter that most affects the results is the search radius for fitting points on a local plane, i.e., the size of the kernel.
To study the effect of this parameter on the results, we segmented the image and evaluated the distresses highlighted when the kernel size changes from 0.1 to 1 m, in 10 cm steps. The stretch of road analyzed is 5 m long (transparent white rectangle in Figure 15). In that stretch there is a pothole (the number 14 visible in Figure 14a) and some shoves and swells (number 44, 45, 46, 47 and 49 in Figure 14b); their geometric parameters are reported in Table 3.
given as input to the software. In particular, in the editing phase, it is important to choose the correct values because the reconstructed 3D model of the road surface represents the reference surface on which the distresses will be highlighted.
The choice of some parameters is based on considerations related to the nature of the input data: the DEM resolution, for example, is in relation to the characteristics of the point cloud, i.e., its density and accuracy. Other parameters are decided according to the type of distress to be identified. For others, some tests have been carried out to verify the influence of the choice of certain values on the results. In the segmentation procedure (see Section 3.3), the parameter that most affects the results is the search radius for fitting points on a local plane, i.e., the size of the kernel.
To study the effect of this parameter on the results, we segmented the image and evaluated the distresses highlighted when the kernel size changes from 0.1 to 1 m, in 10 cm steps. The stretch of road analyzed is 5 m long (transparent white rectangle in Figure 15). In that stretch there is a pothole (the number 14 visible in Figure 14a) and some shoves and swells (number 44, 45, 46, 47 and 49 in Figure 14b); their geometric parameters are reported in Table 3.   Figure 15 shows the results obtained when the kernel size changed. Segmented areas (i.e., those defined as having significant distress) are outlined in black. The color scale indicates the estimated roughness value ∆h (see the legend), which varies between 0 to 2 cm for shoves/swells and 0 to 4 cm for potholes.
Changes in the shape and size of segmented distresses are clearly visible; one can see how distresses tend to appear and increase in size as the kernel size increases. More in details for potholes note that, with a kernel size of 0.10 m, nothing is highlighted because the local plane is too small and follows the trend of the road surface; increasing the kernel size to 0.20 m a region characterized by positive ∆h values (i.e., a pothole) appears; this area stabilizes from 0.60 m onwards, when its size is completely within the kernel size. For shoves/swells, note that, by increasing the kernel size, they gradually appear as they do for potholes but never stabilize.
The results of these tests confirm for us that the choice of value for this parameter is not trivial, hence it is tricky to recommend a criteria of unambiguous choice. A more reliable assessment could be made by analyzing the variation in the severity level of the distress as the kernel size changes, regardless of the variation in the geometric parameters. For this analysis, we have selected some distresses, chosen because they are isolated and well distinguishable, i.e., not characterized by overlapping of different degrades. Figure 16 shows how the level of severity attributed to the analyzed distresses may vary with the size of the kernel. It can be seen that some distress do not change class as the kernel size changes, while others vary more classes as it changes. These last cases concern, in particular, several distresses with positive ∆h value that, more than real potholes, are pitting that smoothly degrade.
for potholes. Changes in the shape and size of segmented distresses are clearly visible; one can see how distresses tend to appear and increase in size as the kernel size increases. More in details for potholes note that, with a kernel size of 0.10 m, nothing is highlighted because the local plane is too small and follows the trend of the road surface; increasing the kernel size to 0.20 m a region characterized by positive h values (i.e., a pothole) appears; this area stabilizes from 0.60 m onwards, when its size is completely within the kernel size. For shoves/swells, note that, by increasing the kernel size, they gradually appear as they do for potholes but never stabilize.
The results of these tests confirm for us that the choice of value for this parameter is not trivial, hence it is tricky to recommend a criteria of unambiguous choice. A more reliable assessment could be made by analyzing the variation in the severity level of the distress as the kernel size changes, regardless of the variation in the geometric parameters. For this analysis, we have selected some distresses, chosen because they are isolated and well distinguishable, i.e., not characterized by overlapping of different degrades. Figure 16 shows how the level of severity attributed to the analyzed distresses may vary with the size of the kernel. It can be seen that some distress do not change class as the kernel size changes, while others vary more classes as it changes. These last cases concern, in particular, several distresses with positive h value that, more than real potholes, are pitting that smoothly degrade. Note also that all swells/shoves analyzed, from a kernel size of 0.60 m onwards, maintain the same level of severity while the potholes vary in severity level, but also tend to stabilize from 0.60 m onwards. Figure 17 graphically shows the values of the most relevant geometric parameters of the distresses listed in Figure 16, as the kernel size changes. In particular, as for the h max value, graphs show that, in most cases, there is a variation in the slope of the lines interpolating the values of h max near the kernel size value of 60 cm. To run the linear interpolation, the 'ransac' function, implemented in Matlab, has been used; the algorithm works, also here, using the MSAC method.
For this reason, in our test area, where the road profile is very irregular, we have chosen a kernel size of 0.60 m. Yet, we want to underline once again the non-triviality of the choice of this parameter. Note also that all swells/shoves analyzed, from a kernel size of 0.60 m onwards, maintain the same level of severity while the potholes vary in severity level, but also tend to stabilize from 0.60 m onwards. Figure 17 graphically shows the values of the most relevant geometric parameters of the distresses listed in Figure 16, as the kernel size changes. In particular, as for the h max value, graphs show that, in most cases, there is a variation in the slope of the lines interpolating the values of h max near the kernel size value of 60 cm. To run the linear interpolation, the 'ransac' function, implemented in Matlab, has been used; the algorithm works, also here, using the MSAC method.
For this reason, in our test area, where the road profile is very irregular, we have chosen a kernel size of 0.60 m. Yet, we want to underline once again the non-triviality of the choice of this parameter.
Another critical aspect about the ∆h computation procedure, is the so-called 'boundary effect', In the roadside areas, for an extension equal to the size of the kernel, the research sphere falls in part out of the roadway, thus the number of points that will be fitted on a plane decrease as you approach the edge of the road (to become about half).
As a result, the ∆h values associated with the points belonging to the boundaries are less reliable than the others. To check this aspect, a simulation was made by computing the ∆h for the points of a central portion of the roadway of about 2 m 2 , selected in the middle of the roadway. Figure 18a shows an excerpt of the point cloud, the grayscale is proportional to the ∆h values computed with a kernel size of 0.60 m. The area framed in white shows the 2 m 2 sub-cloud used for the simulation. For that piece, the difference between the ∆h values computed considering the whole point cloud of Figure 18 (roughness values not influenced by the boundary effect) and the values computed considering only the sub-cloud (values influenced by the boundary effect) has been computed. Another critical aspect about the h computation procedure, is the so-called 'boundary effect', In the roadside areas, for an extension equal to the size of the kernel, the research sphere falls in part out of the roadway, thus the number of points that will be fitted on a plane decrease as you approach the edge of the road (to become about half).
As a result, the h values associated with the points belonging to the boundaries are less reliable than the others. To check this aspect, a simulation was made by computing the Δh for the points of a central portion of the roadway of about 2 m 2 , selected in the middle of the roadway. Figure 18a shows an excerpt of the point cloud, the grayscale is proportional to the h values computed with a kernel size of 0.60 m. The area framed in white shows the 2 m 2 sub-cloud used for the simulation. For that piece, the difference between the h values computed considering the whole point cloud of Figure 18 (roughness values not influenced by the boundary effect) and the values computed considering only the sub-cloud (values influenced by the boundary effect) has been computed. Figure 18b shows the differences h between the two values, classified in 0.5 mm steps (see legend).
As expected, the boundary effect exists, the h in absolute value reaches 3 mm. Hence, the roughness values computed on the boundary strips (outside the area framed in black) are unreliable and there could be a mistake in the attribution of the severity class of the distresses in these zones. It is therefore suggested to exclude from the analysis the whole boundary strip having the width equal to the kernel size.  Figure 18b shows the differences δh between the two values, classified in 0.5 mm steps (see legend).
As expected, the boundary effect exists, the δh in absolute value reaches 3 mm. Hence, the roughness values computed on the boundary strips (outside the area framed in black) are unreliable and there could be a mistake in the attribution of the severity class of the distresses in these zones. It is therefore suggested to exclude from the analysis the whole boundary strip having the width equal to the kernel size.
Remote Sens. 2020, 12, x FOR PEER REVIEW 21 of 24 Figure 18. 'Boundary effect'. (a) A sample of point cloud belonging to a central stretch of road; the area squared in white is the one subject to analysis; the grayscale is proportional to the roughness value h. (b) Classified map of the difference in roughness h computed for the points in the white square of (a). The area within the black box is not influenced by the boundary effect, the area outside the box it is.

Conclusions
In this note, we wanted to highlight the potential of using the LiDAR survey technique, one of the most interesting and continuously developing remote sensing techniques, in order to assess the 'state of degradation' of the road pavement. The technique has proved to be effective in providing data that allow usto build a 3D model of the infrastructure surface.
As the applications of this technique to infrastructure surveying are recent, there are still many critical aspects to evaluate. Commercial software available implement only a few functions, in response to the specific needs of professionals and, therefore, do not fully satisfy the scientific community. In an attempt to overcome these limitations, the problems that arise in the application of the technique to road survey have been analyzed and the studied algorithms have been (b) Classified map of the difference in roughness δh computed for the points in the white square of (a). The area within the black box is not influenced by the boundary effect, the area outside the box it is.

Conclusions
In this note, we wanted to highlight the potential of using the LiDAR survey technique, one of the most interesting and continuously developing remote sensing techniques, in order to assess the 'state of degradation' of the road pavement. The technique has proved to be effective in providing data that allow usto build a 3D model of the infrastructure surface.
As the applications of this technique to infrastructure surveying are recent, there are still many critical aspects to evaluate. Commercial software available implement only a few functions, in response to the specific needs of professionals and, therefore, do not fully satisfy the scientific community. In an attempt to overcome these limitations, the problems that arise in the application of the technique to road survey have been analyzed and the studied algorithms have been implemented in proprietary software (in Matlab ® environment) and tested on a case study.
The point clouds acquired are not directly usable, but they must first be filtered to remove outliers (vegetation, cars, pedestrians, etc.) and extract only the points belonging to the roadway in a semi-automatic way. The implemented algorithm is based on the MSAC method, a robust estimate of the parameters of a mathematical model from a data set containing outliers.
The evaluation of localized distresses in 3D space is rather complex; there is no exact reference surface to estimate any anomalies present with respect to it. In order to quantify and analyze the different types of distress, in particular, the regularity that affects the comfort of users and some localized distresses, mainly responsible for the safety of the infrastructure, it was decided to use a local reference plan, built around each point belonging to the road surface, with respect to which to compute the roughness of each point of the cloud.
The results obtained with the proposed methodology partly depend on the choice of the parameters provided as input to the computation software. The choice of some of them is based on considerations related to the nature of the input data. The resolution of the DEM, for example, is mainly related to the characteristics of the point cloud (density and accuracy) whereas the values to be given to other parameters are suggested by the type of distress to be segmented.
In particular, the choice of kernel size has significant consequences and strongly affects the segmentation process. The sensitivity with which the geometric parameters of the identified regions vary, as the kernel size varies, does not allow an easy and justified choice of the latter. A better evaluation can be made evaluating the variations in the severity levels in place of analyzing the variations in the geometric parameters.
With regard to the procedures implemented, the experimental results have shown that the choice of parameters to be entered by the user for the extraction of the roadway or distress analyzed is not trivial as it strongly influences the results. Due to the complexity of the configuration of the different distresses and the geometry of the roadway, it is appropriate to choose these parameters according to the boundary conditions, it is not possible to identify and set default values valid for any condition.
The segmentation procedure implemented has however allowed us to extract some types of distress in an effective way; the parameters used in the test case have made it possible to correctly segmented and identify the types of distress analyzed, also congruent with some in-situ inspections.