Range image technique for change analysis of rock slopes using dense point cloud data

: The use of a terrestrial laser scanner is examined to measure the changes of rock slopes subject to a wave attack test. Real scenarios are simulated in a water ﬂume facility using a wave attack experiment representing a storm of 3000 waves. The stability of two rock slopes of different steepness was evaluated under the set conditions. For quantiﬁcation of the changes of the slopes after the wave attack test, terrestrial laser scanning was used to acquire dense 3D point cloud data sampling for slope geometries before and after the wave attack experiment. After registration of the two scans, representing situations before and after the wave attack, the cloud-to-cloud distance was determined to identify areas in the slopes that were affected. Then, a range image technique was introduced to generate a raster image to facilitate a change analysis. Using these raster images, volume change was estimated as well. The results indicate that the area around the artiﬁcial coast line is most strongly affected by wave attacks. Another interesting phenomenon considers the change in transport direction of the rocks between the two slopes: from seaward transport for the steeper slope to landward transport for the milder slope. Using the range image technique, the work in this article shows that terrestrial laser scanning is an effective and feasible method for change analysis of long and narrow rock slopes.


Introduction
The protection of coastal zones is fundamental to the effective operation of local socioeconomic activities, especially in the Netherlands, a lowland country with significant coast length. For example, the risks of coastal erosion are expected to increase [1] due to predicted climate change the associated sea level rise, tidal effects, and expected increases in the intensity and frequency of storms. Conversely, coastal erosion poses a serious threat to the nature hazard agents. Various types of coastal elements, i.e., sand banks, sandy beaches, and dunes [2,3], have been used to protect coastal areas from erosion, most of which are made by rockfill [4][5][6][7][8]. With the rapid erosion of coastal structures, demand for natural rocks has increased around the world [9].
A problem that has drawn the attention of researchers involved in coastal protection is which configuration of rocks provides optimal protection in a given wave climate scenario all the time [10]. An extensive review of laboratory studies employing random waves was performed on the stability of photogrammetric point cloud to identify discontinuity sets of rock slopes [49] and the techniques have been found to deliver reliable results.
We discussed an earlier application of TLS for change analysis of rock slopes in [55]. That framework included registration of two point clouds, coordinate transformation from the TLS coordinate system to the so-called slope coordinate system, range image creation, and change detection. However, this framework has two main drawbacks. First, the applied registration and coordinate transformation will introduce small alignment errors, and even minor errors may lead to erroneous changes in the detection results, especially when the magnitude of the change value is relatively small. Second, different pixels in the raster image have identical size but the number of points per pixel has significant differences. In other words, the point density varies which results in local variations in change detection accuracy. To solve these drawbacks related to misalignment and non-uniform accuracy, we propose to introduce the range image technique for change analysis. In this article, the use of a laser scanner to obtain change information is examined and its applicability is confirmed. The aim of the present article is to understand how coastal structures with gentle slopes behave during a wave attack. First, the wave attack test and the scan data acquisition are discussed. Next, to guarantee the accuracy of the change analysis, point cloud registration before and after the wave attack test is performed using the control points method. To identify the area within the entire slope where significant change occurred, the cloud-to-cloud distances between the different slope point clouds are estimated using the CloudCompare software [56]. To obtain further insight in the changes in the affected area, a so-called range image is created based on the vertical and horizontal angular resolutions of the point clouds. Subsequently, quantification of changes is carried out on the two slopes to see how the rock slopes behave against the wave attack.
In this article, first the wave flume and the wave attack tests are described in Section 2.1, and how the dense point cloud is acquired using a Leica C10 scanner. Next, Section 2.2 explains the methodology performed in this article, which is composed of five main steps: (1) point cloud registration between two epochs; (2) damage area detection using cloud-to-cloud distances; (3) range image creation for the interest area; (4) quantification of changes; and (5) volume change estimation. In Section 3, the results are presented and analyzed. Afterwards, a discussion of the effectiveness of the method is provided. Finally, in Section 4 conclusions and future work are presented.

Wave Attack Simulation
In the experiment considered, wave attacks on two types of rock slopes (V:H) were performed, 1:10 and 1:5, respectively. These two types of slopes were assessed considering their safety against wave attack as a guide for the construction of coastal structures. The slopes were tested under set wave conditions to determine their stability in the water flume of Delft University of Technology, which has a length of 40 m, width of 0.8 m, and a height of 0.9 m. Their stability on its turn is derived from the profile change after a wave series that present a Dutch storm event including 3000 waves. The wave heights for slopes of 1:10 and 1:5 are 0.129 m and 0.054 m, respectively. Since the experiment was conducted in a laboratory, artificial lighting from the ceiling was presented. At each occasion, the slope was scanned before the wave attack, and afterwards, after the water was removed from the flume. To mount the terrestrial laser scanner directly above the rock slope, a special wooden construction was made, see Figure 1a. In the tests, the placed wooden construction was regarded stable during each scan. To simulate real coastal conditions, the wave surface did not cover the entire slope as can be seen in Figure 2, where oblique surfaces surrounded by solid lines are rock "Slope 1" (1:10) and rock "Slope 2" (1:5), respectively. The blue dashed line represents the wave surface in the still state. The red dashed lines indicate the intersection line (coast line) between still wave level and the two rock slopes, respectively. The water flume was designed for modeling coastal structures and assessing their safety against wave attack, storm surges, and flooding. As shown in Figure 2, a special wave maker was set at the right side of the water flume. The wave generator paddle at a 0° angle wave attack provided mechanically irregular waves on the flume. The generator was controlled by a computer program [57]. The program was developed for simulating regular, irregular, and a variety of other wave conditions. In our experiment, a Dutch storm of 3000 waves was generated by the JONWAP (Joint North Sea Wave Project) wave spectrum.
In this experiment, two mild slopes (1:10 and 1:5) were tested under set wave conditions to determine their stability. The stones used on the slopes have a nominal diameter of 16.20 mm and were placed at a layer thickness of two times the nominal diameter (i.e., 32.40 mm). The shape of the stone were determined by their length-to-thickness ratio, their blockiness, their roundness, or their cubocity, as described in the Rock Manual [58]. To simulate real coastal conditions, the wave surface did not cover the entire slope as can be seen in Figure 2, where oblique surfaces surrounded by solid lines are rock "Slope 1" (1:10) and rock "Slope 2" (1:5), respectively. The blue dashed line represents the wave surface in the still state. The red dashed lines indicate the intersection line (coast line) between still wave level and the two rock slopes, respectively. To simulate real coastal conditions, the wave surface did not cover the entire slope as can be seen in Figure 2, where oblique surfaces surrounded by solid lines are rock "Slope 1" (1:10) and rock "Slope 2" (1:5), respectively. The blue dashed line represents the wave surface in the still state. The red dashed lines indicate the intersection line (coast line) between still wave level and the two rock slopes, respectively. The water flume was designed for modeling coastal structures and assessing their safety against wave attack, storm surges, and flooding. As shown in Figure 2, a special wave maker was set at the right side of the water flume. The wave generator paddle at a 0° angle wave attack provided mechanically irregular waves on the flume. The generator was controlled by a computer program [57]. The program was developed for simulating regular, irregular, and a variety of other wave conditions. In our experiment, a Dutch storm of 3000 waves was generated by the JONWAP (Joint North Sea Wave Project) wave spectrum.
In this experiment, two mild slopes (1:10 and 1:5) were tested under set wave conditions to determine their stability. The stones used on the slopes have a nominal diameter of 16.20 mm and were placed at a layer thickness of two times the nominal diameter (i.e., 32.40 mm). The shape of the stone were determined by their length-to-thickness ratio, their blockiness, their roundness, or their cubocity, as described in the Rock Manual [58]. The water flume was designed for modeling coastal structures and assessing their safety against wave attack, storm surges, and flooding. As shown in Figure 2, a special wave maker was set at the right side of the water flume. The wave generator paddle at a 0 • angle wave attack provided mechanically irregular waves on the flume. The generator was controlled by a computer program [57]. The program was developed for simulating regular, irregular, and a variety of other wave conditions. In our experiment, a Dutch storm of 3000 waves was generated by the JONWAP (Joint North Sea Wave Project) wave spectrum.
In this experiment, two mild slopes (1:10 and 1:5) were tested under set wave conditions to determine their stability. The stones used on the slopes have a nominal diameter of 16.20 mm and were placed at a layer thickness of two times the nominal diameter (i.e., 32.40 mm). The shape of the stone were determined by their length-to-thickness ratio, their blockiness, their roundness, or their cubocity, as described in the Rock Manual [58].

Point Cloud Acquisition
The experiment was performed from 14 August 2016 to 28 September 2016. Repeated scanning was done before and after the wave attack testing. A C10 Scan Station scanner (Leica, Heerbrugg, Switzerland) was used for this experiment. The Leica C10 Scan Station scanner, with an effective range of 280 m at 90% reflectively, exploits a time-of-flight principle for measuring the range between the scanner and the object. The accuracy of a single measurement is 6 mm (one sigma) in position and 4 mm (one sigma) in depth at ranges up to 50 m [59]. For each epoch, two scans were made. The first scan was made based on a minimum resolution corresponding to 0.2 m in horizontal and vertical spacing when the range is 100 m [40], and with a complete field of view, i.e., 360 arc-degree with respect to the Z-axis and 270 arc-degree for vertical amplitude. While the second scan was made based on a high-resolution with an area only including the water flume, in which the high-resolution corresponds to 0.05 m in horizontal and vertical spacing at a range of 100 m. Table 1 describes the point clouds acquired for two slopes. The location of the scan position was set on the top side of the slope. To guarantee the accuracy of the registration, three targets with relatively strong geometry were positioned on the water flume frame, see the right figure of Figure 1. In this figure, the positions of the scanner, the targets and the slope are indicated in a scanner view of the experiment scene.

Methodology
The methodology is based on the dense point cloud of the rock slope before and after the wave attack test. A workflow summarizing the methodology is presented in Figure 3 and the steps are described in the following paragraphs. The experiment was performed from 14 August 2016 to 28 September 2016. Repeated scanning was done before and after the wave attack testing. A C10 Scan Station scanner (Leica, Heerbrugg, Switzerland) was used for this experiment. The Leica C10 Scan Station scanner, with an effective range of 280 m at 90% reflectively, exploits a time-of-flight principle for measuring the range between the scanner and the object. The accuracy of a single measurement is 6 mm (one sigma) in position and 4 mm (one sigma) in depth at ranges up to 50 m [59]. For each epoch, two scans were made. The first scan was made based on a minimum resolution corresponding to 0.2 m in horizontal and vertical spacing when the range is 100 m [40], and with a complete field of view, i.e., 360 arc-degree with respect to the Z-axis and 270 arc-degree for vertical amplitude. While the second scan was made based on a high-resolution with an area only including the water flume, in which the high-resolution corresponds to 0.05 m in horizontal and vertical spacing at a range of 100 m. Table 1 describes the point clouds acquired for two slopes. The location of the scan position was set on the top side of the slope. To guarantee the accuracy of the registration, three targets with relatively strong geometry were positioned on the water flume frame, see the right figure of Figure 1. In this figure, the positions of the scanner, the targets and the slope are indicated in a scanner view of the experiment scene.

Methodology
The methodology is based on the dense point cloud of the rock slope before and after the wave attack test. A workflow summarizing the methodology is presented in Figure 3 and the steps are described in the following paragraphs.

Registration of Two Epochs
In general, the locations of the scanning stations were not strictly the same at the different epochs, so the coordinates acquired in two consecutive epochs are not expected to be at the same reference system with different origins and/or orientation. Thus, the registration of point clouds from two epochs is required before change analysis can be performed. Registration aligns and combines multiple point clouds into a single set of range data. However, the introduced registration error will inevitably affect the change detection results, even though the magnitude is relatively small.
In our research, the control points method (TCP) is employed for the process of registration of various scans. Three targets made by special materials manufactured by the Leica Company

Registration of Two Epochs
In general, the locations of the scanning stations were not strictly the same at the different epochs, so the coordinates acquired in two consecutive epochs are not expected to be at the same reference system with different origins and/or orientation. Thus, the registration of point clouds from two epochs is required before change analysis can be performed. Registration aligns and combines multiple point clouds into a single set of range data. However, the introduced registration error will inevitably affect the change detection results, even though the magnitude is relatively small.
In our research, the control points method (TCP) is employed for the process of registration of various scans. Three targets made by special materials manufactured by the Leica Company (Wetzlar, Germany) are used, the centroid of which can be extracted by the scanner itself in a highly accurate way. In this way, the Leica Cyclone software could provide an automatic registration based on TCP. Considering the characteristic of the water flume with long length compared to its width, it is difficult to position the targets in a geometrically strong way. However, due to the high accuracy of target center extraction algorithm by Leica Cyclone software, a high-precision registration could be performed using the Least Squares method: the errors for the three targets between two epochs are less than or equal to 1 mm.

Cloud-to-Cloud Distances
As a clean and proper model of a subject surface is difficult to obtain, the classical way to detect changes from point clouds is cloud-to-cloud distances based on local neighbors [60]. This is a fast and direct 3D comparison method that has been widely applied in point cloud comparison. The method has the advantage that it does not require meshing of the dataset, nor estimation of point normals. For each point in the evaluated point cloud, the nearest point in the reference cloud is searched and the Euclidean distance between them is computed. Direct point to point comparison has at least two disadvantages: results are often noisy, and results are dominated by local point spacing. Alternatively, from one point cloud a surface approximation is estimated. The process is first to find the nearest point in the reference cloud. Based on the nearest point and its predefined number of neighbors, the surface of the reference is modeled locally by fitting a mathematical primitive like a plane. The distance to this local model is finally determined. A graphical diagram illustrating the local of a local surface model for cloud-to-cloud distances is shown in Figure 4. (Wetzlar, Germany) are used, the centroid of which can be extracted by the scanner itself in a highly accurate way. In this way, the Leica Cyclone software could provide an automatic registration based on TCP. Considering the characteristic of the water flume with long length compared to its width, it is difficult to position the targets in a geometrically strong way. However, due to the high accuracy of target center extraction algorithm by Leica Cyclone software, a high-precision registration could be performed using the Least Squares method: the errors for the three targets between two epochs are less than or equal to 1 mm.

Cloud-to-Cloud Distances
As a clean and proper model of a subject surface is difficult to obtain, the classical way to detect changes from point clouds is cloud-to-cloud distances based on local neighbors [60]. This is a fast and direct 3D comparison method that has been widely applied in point cloud comparison. The method has the advantage that it does not require meshing of the dataset, nor estimation of point normals. For each point in the evaluated point cloud, the nearest point in the reference cloud is searched and the Euclidean distance between them is computed. Direct point to point comparison has at least two disadvantages: results are often noisy, and results are dominated by local point spacing. Alternatively, from one point cloud a surface approximation is estimated. The process is first to find the nearest point in the reference cloud. Based on the nearest point and its predefined number of neighbors, the surface of the reference is modeled locally by fitting a mathematical primitive like a plane. The distance to this local model is finally determined. A graphical diagram illustrating the local of a local surface model for cloud-to-cloud distances is shown in Figure 4. One advantage of using cloud-to-cloud distances for change detection is that it does not require the explicit pairing of point in two datasets [61]. However, cloud-to-cloud distances are always positive which means that the direction of change is still unknown.

Range Image and Raster Image Creation
Actually, the scanner mechanism can be considered to operate in a spherical coordinate system. In most current systems, the point cloud acquired is expressed in Cartesian coordinates, which means that the original acquisition organization is lost. Therefore, this unorganized cartesian representation is converted to spherical coordinates before further processing. As a result, the unorganized 3D data are expressed in a 2D range image format. The angular resolution θ Δ is predefined depending on the scanning resolution that can be regarded as inherent parameter of the scan. For searching in a simplified way, a conversion is made to generate regular raster (2D coordinates) [62]:  One advantage of using cloud-to-cloud distances for change detection is that it does not require the explicit pairing of point in two datasets [61]. However, cloud-to-cloud distances are always positive which means that the direction of change is still unknown.

Range Image and Raster Image Creation
Actually, the scanner mechanism can be considered to operate in a spherical coordinate system. In most current systems, the point cloud acquired is expressed in Cartesian coordinates, which means that the original acquisition organization is lost. Therefore, this unorganized cartesian representation is converted to spherical coordinates before further processing. As a result, the unorganized 3D data are expressed in a 2D range image format. The angular resolution ∆θ is predefined depending on the scanning resolution that can be regarded as inherent parameter of the scan. For searching in a simplified way, a conversion is made to generate regular raster (2D coordinates) [62]: where θ and ϕ are the 2D spherical coordinates of a point derived from the 3D coordinates; X, Y, and Z are the 3D Cartesian coordinates as gathered by the scanner; ∆θ H and ∆θ V are the horizontal and vertical angular resolutions, respectively, which depend on the scanner settings. The range from the point to the scanner center is determined by ρ = In such a coordinate system, the 2D coordinates of all points are organized in a 2D array indexed by row and column number.
In contrast, for a point P i (θ i , ϕ i , ρ i ) in the range image, the corresponding Cartesian coordinates can be expressed as Since there is no absolute corresponding point between two epochs of the same scene, the raster image is generated to realize the change in value acquisition by comparing the corresponding grid at the same location in different epochs. Once the point cloud is in this 2D space, the next task consists of converting this dataset, composed of isolated points, into a raster image where the size of the grid in the raster image is predefined by the user, depending the angular resolution and practical necessity. In this raster image, the attribute of each grid is the mean range obtained by averaging all ranges of points in the current gridded. Assuming that the pixel size is k times the angular resolution, the size of the bounding box is k∆θ H × k∆θ V , which denotes the resolution of the raster grid cells. It should be noted that the resolution here is not the real dimensions in 3D space, but the dimensions related to angular resolution in range image as introduced above. In practice, the user can fix the resolution of the pixel size or estimate it depending on the point cloud density in the 2D space. Consequently, a raster grid of size m × n was built, where m is the number of rows and n is the number of columns.

Quantification of Changes
Once the raster images for both point clouds are created, the range change is subsequently tracked by comparing the corresponding cells in two different epochs, which is computed as where ρ 1 and ρ 2 denote the ranges to the scanner center for the first and second epochs, ∆ρ is the range change between two epochs. Although the range change for each cell was estimated, the change along or perpendicular to the slope is still unclear. For further exploiting the change information, it is better to project the range change to the directions along the slope and perpendicular to the slope, the details of this process are as follows.
(1) Slope direction determination. As shown in Figure 1, the scanner is mounted at the middle along the cross section of the trunk. In such a case, the cells at the middle of the water flume along the slope direction have the same horizontal angle in the range image. Subsequently, the horizontal angle (θ M ) representing the slope direction is obtained by averaging the horizontal angles of all cells in the range image. By connecting the cells with the same horizontal angle θ M the slope direction is determined. Afterwards, the 3D points in the current cells representing the slope direction are extracted. Using the points in 3D space, the vectors for the direction are estimated by means of Principle Component Analysis (PCA), the principle of which is shown in Jolliffe [63] in detail. Here, we define the direction away from the scanner as the position direction (down-slope). A vector V 1 is used to represent the slope direction. (2) Direction across the water flume determination. Once the slope direction is determined, the next step is to determine the section direction. It is apparent that the section direction is perpendicular to the slope direction, so the key is to determine the position direction of the section. In our research, the direction from the left to the right from the viewpoint of the scanner is regarded as the positive direction of the section. A vector V 2 is used to represent the section direction. (3) Direction vertically upwards against the slope plane. The direction (V 3 ) is pointing upwards perpendicular to the slope plane determined by the vectors V 1 and V 2 , which is simply given as where '×' denotes cross multiplication.
Once we have defined the changes in the three directions V 1 , V 2 , and V 3 , we could project the range change to the three directions, as presented subsequently.
Indeed, the range for a cell in the range image is the distance from the location of the current cell to the scanner center. In this sense, the range also could be expressed by a vector connecting the current cell to the scanner center. Let P i (X i , Y i , Z i ) denote the 3D Cartesian coordinates for the current cell. As the 3D Cartesian coordinates of the scanner center are always O(0, 0, 0), the range direction is therefore expressed as Subsequently, the range change projection is performed according to Equation (5).
where ∆ρ i is a scalar denoting the range change for the current cell, P i O, V 1 , V 2 , and V 3 are vectors, and • is the norm of a vector.

Volume Change Estimation
Once the range change is obtained, volume change detection is performed based on the cells in the raster image. Given a cell C(i, j)|i = 1, 2, · · · , m; j = 1, 2, · · · , n in the raster image, the 2D coordinates of four corners are determined, i.e., A 1 (θ 1 , ϕ 1 ), A 2 (θ 2 , ϕ 2 ), A 3 (θ 3 , ϕ 3 ), and A 4 (θ 4 , ϕ 4 ). Next, the 3D coordinates of the four corners ( B i (X i , Y i , Z i )|i = 1, 2, 3, 4 ) in Cartesian coordinates are computed according to Equation (1). The cell in 3D space is regarded as a rectangle, the area (S ij )is therefore calculated using the 3D coordinates of four corners. Finally, combined with the change value vertically upwards against the slope plane provided by the previous step, the volume change is estimated as ∆V ij = L z · S ij . Here we should remark that the cell projecting back to the 3D space is not rigorously a rectangle but an annulus sector ( Figure 5), which can be approximated by a trapezoid since the angles ∆θ are close to zero. Namely when ∆θ → 0 , the chord and arc with respect to the identical angle are with nearly equal values of length. As illustrated in Figure 5, the current red cell in raster image is described by four corners A 1 , A 2 , A 3 , and A 4 while the corresponding area in the 3D space is represented by B 1 , B 2 , B 3 , and B 4 , such that In fact, the connection of points B 1 and B 2 is an arc rather than a chord (straight line). The same holds for the connection between points B 3 and B 4 . In addition, the length from B 1 to B 2 is not strictly equal to the length from B 3 to B 4 as they have the same horizontal angular difference but a different vertical angular difference in 3D space. However, the cell size is relatively small, that is, the length difference induced by the vertical angular difference is small enough to be ignored. Given that in our Remote Sens. 2018, 10, 1792 9 of 25 experiment, the maximum range in the scene is less than 20 m, the cell in 3D space is treated as a rectangle in the later process without compromising the final accuracy. 2 4 In fact, the connection of points 1 B and 2 B is an arc rather than a chord (straight line). The same holds for the connection between points 3 B and 4 B . In addition, the length from 1 B to 2 B is not strictly equal to the length from 3 B to 4 B as they have the same horizontal angular difference but a different vertical angular difference in 3D space. However, the cell size is relatively small, that is, the length difference induced by the vertical angular difference is small enough to be ignored. Given that in our experiment, the maximum range in the scene is less than 20 m, the cell in 3D space is treated as a rectangle in the later process without compromising the final accuracy. Figure 5. 3D profiles for scanning, in which the transformation between raster image and 3D coordinates is illustrated: For example, assuming that the range for the current cell to the scanner center is 20 m, the vertical and horizontal angle are 1.00 rad and 1.00 rad, respectively, the horizontal and vertical angular resolutions are 0.0005 rad, and the pixel size in 2D space is set as 0.0025 rad × 0.0025 rad, the 3D coordinates for the four corners in 3D space are therefore ( ) 1

Registration
As mentioned in Section 2.1, three plane targets were set for aligning the point clouds from different epochs. Control points were additionally measured by the scanner, to align the scans more accurately to a common reference system for the coordinates. The errors for the registration in Slope 1 are 0.000, 0.001, and 0.001 mm along the X, Y, and Z axes, respectively. While in Slope 2, the errors are 0.001, 0.000, and 0.000 mm along the X, Y, and Z axes, respectively. These errors indicate that the registration was performed with high precision for two slopes, which therefore guarantees the reliability of change detection results in the same reference system.

Change Analysis Using Cloud-to-Cloud Distances
Once two 3D point clouds are spatially registered for two slopes, they can be compared against each other for change analysis. Thus, the cloud-to-cloud distances for the two slopes were . 3D profiles for scanning, in which the transformation between raster image and 3D coordinates is illustrated: For example, assuming that the range for the current cell to the scanner center is 20 m, the vertical and horizontal angle are 1.00 rad and 1.00 rad, respectively, the horizontal and vertical angular resolutions are 0.0005 rad, and the pixel size in 2D space is set as 0.0025 rad × 0.0025 rad, the 3D coordinates for the four corners in 3D space are therefore B 1 (9.093, 5.8385, 16.8294), B 2 (9.1075, 5.8158, 16.8294), B 3 (9.0575, 5.8158, 16.8564), and B 4 (9.0721, 5.7931, 16.8564). The difference between B 1 B 2 and B 3 B 4 is then estimated, which is equal to 0.0001 m.

Registration
As mentioned in Section 2.1, three plane targets were set for aligning the point clouds from different epochs. Control points were additionally measured by the scanner, to align the scans more accurately to a common reference system for the coordinates. The errors for the registration in Slope 1 are 0.000, 0.001, and 0.001 mm along the X, Y, and Z axes, respectively. While in Slope 2, the errors are 0.001, 0.000, and 0.000 mm along the X, Y, and Z axes, respectively. These errors indicate that the registration was performed with high precision for two slopes, which therefore guarantees the reliability of change detection results in the same reference system.

Change Analysis Using Cloud-to-Cloud Distances
Once two 3D point clouds are spatially registered for two slopes, they can be compared against each other for change analysis. Thus, the cloud-to-cloud distances for the two slopes were calculated using the CloudCompare software [56] and the resulting scatter diagrams colored by distance value are shown in Figure 6 (for slope 1) and Figure 7 (for slope 2), respectively. As the stones on the slope do not have regular shapes, a least squares best fitting plane may not always be a good approximation. Therefore, quadric fitting of local surfaces is employed in our study. When searching for the k-nearest neighbors (k-NN) points, six nearby points were used. In the case of Slope 1 (1:10), the maximum cloud-to-cloud distance was 0.0333 m, while for Slope 2 (1:5) the value was 0.0206 m. The cloud to cloud distances variation distribution is summarized for Slope 1 and Slope 2 respectively, see Tables 2 and 3. the slope do not have regular shapes, a least squares best fitting plane may not always be a good approximation. Therefore, quadric fitting of local surfaces is employed in our study. When searching for the k-nearest neighbors (k-NN) points, six nearby points were used. In the case of Slope 1 (1:10), the maximum cloud-to-cloud distance was 0.0333 m, while for Slope 2 (1:5) the value was 0.0206 m. The cloud to cloud distances variation distribution is summarized for Slope 1 and Slope 2 respectively, see Tables 2 and 3.     approximation. Therefore, quadric fitting of local surfaces is employed in our study. When searching for the k-nearest neighbors (k-NN) points, six nearby points were used. In the case of Slope 1 (1:10), the maximum cloud-to-cloud distance was 0.0333 m, while for Slope 2 (1:5) the value was 0.0206 m. The cloud to cloud distances variation distribution is summarized for Slope 1 and Slope 2 respectively, see Tables 2 and 3.       For Slope 1, the number of points with a cloud-to-cloud distance less than 0.003 mm is 2,505,973, which accounts for 93.7% of all points of the slope (Table 2). While in the case of Slope 2, this percentage is 97.8% as can be inferred from Table 3. This indicates that a relatively small change occurred at most of the surface of each slope. As seen from Figures 6 and 7, larger changes occurred near the intersection line for both slopes, which indicates the region most affected by the wave attack. A possible explanation for this phenomenon is that the location near the intersection line is influenced significantly by the series of waves, and the runup and rundown velocities in this location are highest compared to the other areas. By estimating the cloud-to-cloud distances, the affected area and the associated change values could be identified. Noteworthy, the cloud-to-cloud distances are absolute values, and therefore the direction of change cannot be detected. That is, in our case the erosion or dilation for the rock slopes cannot be distinguished.

Range Image Technique for Quantifying Change Regions
From change analysis results using the cloud-to-cloud distance, it is apparent that the change in the area near the intersection line is larger than in other areas. Therefore, points were selected manually from the global point cloud around this area. The range images were created according to Equation (1) for the two slopes respectively. Figure 8a,c show the raw point cloud colored by cloud-to-cloud distances for Slope 1 and Slope 2, respectively. Figure 8b,d show the range images colored by ranges for Slope 1 and Slope 2, respectively. For Slope 1, the number of points with a cloud-to-cloud distance less than 0.003 mm is 2,505,973, which accounts for 93.7% of all points of the slope ( Table 2). While in the case of Slope 2, this percentage is 97.8% as can be inferred from Table 3. This indicates that a relatively small change occurred at most of the surface of each slope. As seen from Figures 6 and 7, larger changes occurred near the intersection line for both slopes, which indicates the region most affected by the wave attack. A possible explanation for this phenomenon is that the location near the intersection line is influenced significantly by the series of waves, and the runup and rundown velocities in this location are highest compared to the other areas. By estimating the cloud-to-cloud distances, the affected area and the associated change values could be identified. Noteworthy, the cloud-to-cloud distances are absolute values, and therefore the direction of change cannot be detected. That is, in our case the erosion or dilation for the rock slopes cannot be distinguished.

Range Image Technique for Quantifying Change Regions
From change analysis results using the cloud-to-cloud distance, it is apparent that the change in the area near the intersection line is larger than in other areas. Therefore, points were selected manually from the global point cloud around this area. The range images were created according to Equation (1) for the two slopes respectively. Figure 8a,c show the raw point cloud colored by cloud-to-cloud distances for Slope 1 and Slope 2, respectively. Figure 8b,d show the range images colored by ranges for Slope 1 and Slope 2, respectively.  In the range image, the horizontal and vertical axes correspond to the horizontal and vertical scan angles, respectively. In our experiment, the horizontal and vertical angle resolutions are both 0.0005 rad. Here, the size of the image cell was set at 5 times this angular resolution, i.e., at a value of 0.0025 rad, to insure the number of points in each cell is more than 20. Next, the range image was divided into cells of 0.0025 rad × 0.0025 rad. To remove the influence of edge effects, cell creation started from the center of the range image. For each point in a given cell, the mean range was estimated by averaging the ranges of all point within the current cell. Through comparing the corresponding cells in two epochs, the range change was estimated. The range images colored by range change for Slope 1 and Slope 2 are shown in Figures 9 and 10, respectively. Range changes with negative values indicate rocks in the current cell that have moved closer to the scanner location after the wave attack test. Conversely, the range change with positive values correspond to the rocks in the current cell that moved away from the scanner center after the wave attack test. of 0.0025 rad, to insure the number of points in each cell is more than 20. Next, the range image was divided into cells of 0.0025 rad × 0.0025 rad. To remove the influence of edge effects, cell creation started from the center of the range image. For each point in a given cell, the mean range was estimated by averaging the ranges of all point within the current cell. Through comparing the corresponding cells in two epochs, the range change was estimated. The range images colored by range change for Slope 1 and Slope 2 are shown in Figures 9 and 10, respectively. Range changes with negative values indicate rocks in the current cell that have moved closer to the scanner location after the wave attack test. Conversely, the range change with positive values correspond to the rocks in the current cell that moved away from the scanner center after the wave attack test.  For further characterizing the change information, the range change is projected to the directions along and perpendicular to the slope. Based on the horizontal angles, the center of the slope is obtained according to the method introduced in Section 3.3. Afterwards, the vectors representing the main direction ( 1 V ) along the slope, the second direction ( 2 V ) across the water flume, and the third direction ( 3 V ) vertically upwards against the slope plane are estimated. These vectors are summarized in Table 4. divided into cells of 0.0025 rad × 0.0025 rad. To remove the influence of edge effects, cell creation started from the center of the range image. For each point in a given cell, the mean range was estimated by averaging the ranges of all point within the current cell. Through comparing the corresponding cells in two epochs, the range change was estimated. The range images colored by range change for Slope 1 and Slope 2 are shown in Figures 9 and 10, respectively. Range changes with negative values indicate rocks in the current cell that have moved closer to the scanner location after the wave attack test. Conversely, the range change with positive values correspond to the rocks in the current cell that moved away from the scanner center after the wave attack test.  For further characterizing the change information, the range change is projected to the directions along and perpendicular to the slope. Based on the horizontal angles, the center of the slope is obtained according to the method introduced in Section 3.3. Afterwards, the vectors representing the main direction ( 1 V ) along the slope, the second direction ( 2 V ) across the water flume, and the third direction ( 3 V ) vertically upwards against the slope plane are estimated. These vectors are summarized in Table 4. For further characterizing the change information, the range change is projected to the directions along and perpendicular to the slope. Based on the horizontal angles, the center of the slope is obtained according to the method introduced in Section 3.3. Afterwards, the vectors representing the main direction (V 1 ) along the slope, the second direction (V 2 ) across the water flume, and the third direction (V 3 ) vertically upwards against the slope plane are estimated. These vectors are summarized in Table 4. Subsequently, the range change for every cell in the range image is decomposed in the three directions as discussed above. The histograms of the number of cells along slope, across slope, and upwards for Slope 1 are shown in Figure 11a-c, respectively. Similarly, the histograms of the number of cells along slope, across slope, and upwards for Slope 2 are shown in Figure 11d-f, respectively. Figure 11. Histogram of changes. Panels (a-c) are the histograms for Slope 1 in the directions along the slope, across the water flume direction, and vertically upwards against the slope plane, respectively. Panels (d-f) are the histograms for Slope 2 in the directions along the slope, across the water flume direction, and vertically upwards against the slope plane, respectively.
As can be seen in Figure 11 In addition, the cells with an absolute change value less than 0.001 m in the direction across the water flume occupy 92.5% and 95.7% of Slope 1 and Slope 2, respectively. In this regard, considering the registration error, the surface density in different scans, etc., the change in the direction across the water flume can be ignored in the further analysis. Next, the scatter diagrams for the two slopes in the direction along slope are plotted, in Figures 12 and 13.
As shown in Figures 12 and 13, the changes in the direction along slope are symmetrically distributed along the center line of the slope length, i.e., the center line of the water flume. Thus, the change profile along the slope length, i.e., the vertical axis direction (ϕ) in 2D space, is defined and computed by averaging all the cell values within the same horizontal axis (θ), as in Equation (6): where n is the number of columns, in the raster image, and C(i, j) denotes the change of a raster grid in the i-th row and j-th column.   As shown in Figures 12 and 13, the changes in the direction along slope are symmetrically distributed along the center line of the slope length, i.e., the center line of the water flume. Thus, the change profile along the slope length, i.e., the vertical axis direction ( ϕ ) in 2D space, is defined and computed by averaging all the cell values within the same horizontal axis (θ ), as in Equation (6)   As shown in Figures 12 and 13, the changes in the direction along slope are symmetrically distributed along the center line of the slope length, i.e., the center line of the water flume. Thus, the change profile along the slope length, i.e., the vertical axis direction ( ϕ ) in 2D space, is defined and computed by averaging all the cell values within the same horizontal axis (θ ), as in Equation (6) Table 5.
for Slope 2 are 0.0052 m and −0.0032 m, respectively, while the same change values for Slope 1 are more significant with values of 0.0121 m and −0.0132 m, respectively. Compared to the change in the direction along the slope, the changes in the direction vertically upwards against the slope plane are less significant, i.e., the maximum change values for Slope 1 and Slope 2 are 0.0036 m and 0.001 m, respectively, and the minimum change values are −0.0036 m and −0.0018 m, respectively. The change information is also summarized in Table 5.  Close examination of the change information in Figure 14 and Table 5 shows that most change happened at the area around the intersection line. In addition, the changes in the direction along the slope have a similar tendency for Slope 1 and Slope 2, i.e., along the slope, first the change becomes  Panels (a,b) show the change in the directions along the slope length and vertically upwards against the slope plane for slope 1, while panels (c,d) show the change in the directions along the slope length and vertically upwards against the slope plane respectively for slope 2. Each value in the profile is obtained by averaging the changes (i.e., along the slope length or across the water flume) of cells in the same rows. Close examination of the change information in Figure 14 and Table 5 shows that most change happened at the area around the intersection line. In addition, the changes in the direction along the slope have a similar tendency for Slope 1 and Slope 2, i.e., along the slope, first the change becomes larger, which indicates that some materials (stones) moving to this location from other areas. Then the change becomes smaller which indicates the materials here are moving away after the wave attack test. The slopes tend to be stable away from this location. However, the changes vertically upwards against the slope plane act differently for two slopes. That is, for Slope 1, the change is first larger, then smaller. For Slope 2, the change first is smaller, then larger. This is an interesting phenomenon and the following conclusions can be drawn, (1) For the steeper 1:5 slope, the transport direction changes to be mostly seaward, that is, towards the direction away from the scanner and (2) for the less steep 1:10 slope, the transport is in the opposite direction, that is, towards the direction close to the scanner.
In our previous work [55], we analyzed the same datasets for the slopes but performed another approach consisting of exploiting a slope coordinate system, coordinate transformation, range image creation, change detection, etc. Through comparing the results, we found that the results acquired by the approach proposed in this article correspond well to the results presented in our previous contribution [55] including the change values, change tendency, etc.

Volume Change Estimation Results
After obtaining the range change for the area of interest, volume change detection is performed according to the method introduced in Section 2.2.4. Figure 15a,b show the histogram and scatter diagram of volume change for Slope 1 while Figure 16a,b show the histogram and scatter diagram of volume change for Slope 2. The maximum volume change, i.e., erosion that occurred for Slope 1 and Slope 2 are 1257.1 mm 3 and 1004.8 mm 3 , respectively. The minimum volume change, i.e., dilation that occurred for Slope 1 and Slope 2 are 1229.8 mm 3 and 528.5 mm 3 , respectively. By summing all the cells in the analyzed area, the total volume changes are estimated with the values of −79,060.9 mm 3 and −22,552.2 mm 3 for Slope 1 and Slope 2, respectively. It is evident that the volume change in Slope 1 is greater than the volume change in Slope 2. This partly reveals that the influence induced by the wave attack to Slope 1 is greater than to Slope 2, which caused by the different wave conditions, i.e., the wave height applied to slope 1 (0.129 m) was larger than the one that was used for slope 2 (0.054 m).
phenomenon and the following conclusions can be drawn, (1) For the steeper 1:5 slope, the transport direction changes to be mostly seaward, that is, towards the direction away from the scanner and (2) for the less steep 1:10 slope, the transport is in the opposite direction, that is, towards the direction close to the scanner.
In our previous work [55], we analyzed the same datasets for the slopes but performed another approach consisting of exploiting a slope coordinate system, coordinate transformation, range image creation, change detection, etc. Through comparing the results, we found that the results acquired by the approach proposed in this article correspond well to the results presented in our previous contribution [55] including the change values, change tendency, etc.

Volume Change Estimation Results
After obtaining the range change for the area of interest, volume change detection is performed according to the method introduced in Section 2.2.4. Figure 15a,b show the histogram and scatter diagram of volume change for Slope 1 while Figure 16a,b show the histogram and scatter diagram of volume change for Slope 2. The maximum volume change, i.e., erosion that occurred for Slope 1 and Slope 2 are 1257.1 mm 3 and 1004.8 mm 3 , respectively. The minimum volume change, i.e., dilation that occurred for Slope 1 and Slope 2 are 1229.8 mm 3 and 528.5 mm 3 , respectively. By summing all the cells in the analyzed area, the total volume changes are estimated with the values of −79,060.9 mm 3 and −22,552.2 mm 3 for Slope 1 and Slope 2, respectively. It is evident that the volume change in Slope 1 is greater than the volume change in Slope 2. This partly reveals that the influence induced by the wave attack to Slope 1 is greater than to Slope 2, which caused by the different wave conditions, i.e., the wave height applied to slope 1 (0.129 m) was larger than the one that was used for slope 2 (0.054 m).

Limitation Analysis
In this study, the rock slopes considered are long and narrow which ensures the excellent performance of the range image technique. However, the width of the slope influences the change profile along the slope length generation. The typical graphical illustration of the scanning is described in Figure 17 where O denotes the scanner center and O′ denotes the scanner center projection in the slope plane. The slope length direction is represented by O C ′ , the points in the arc ACB are with identical vertical angle so they are expected to have the same vertical coordinates in 2D space that are revealed by a line. In reality, the points in the line AB across the slope are expected to be extracted rather than the point in the arc ACB . This somehow causes errors in generating the change profile along the slope length. Assuming that the distance from O′ to D in the slope plane is ′ with a value of L , the distance from the scanner center to the

Limitation Analysis
In this study, the rock slopes considered are long and narrow which ensures the excellent performance of the range image technique. However, the width of the slope influences the change profile along the slope length generation. The typical graphical illustration of the scanning is described in Figure 17 where O denotes the scanner center and O denotes the scanner center projection in the slope plane. The slope length direction is represented by O C, the points in the arc ACB are with identical vertical angle so they are expected to have the same vertical coordinates in 2D space that are revealed by a line. In reality, the points in the line AB across the slope are expected to be extracted rather than the point in the arc ACB. This somehow causes errors in generating the change profile along the slope length. Assuming that the distance from O to D in the slope plane is O D with a value of L 1 , the distance from the scanner center to the slope plane is OO with a value of H, the vertical angle at current position is ϕ, the angle between OA and OD is θ, and the width of the slope (AB) is W 1 , thus, In our experiment, the width of the slope is 0.8 m, the range for Slope 1 is less than 9 m while for Slope 2 it is less than 15 m. The distance from the scanner center to the slope plane is less than 2 m. Here we define the distance OO′ equals to 2 m, and the range of 2 m to 15 m is used to illustrate the influence induced by the range image technique. A graph indicating the difference ( C D ) against the ranges is shown in Figure 18. It is apparently that the difference decreases as the range increase. The difference at the range of 2 m is 0.0396 m while the difference at the range of 15 m is 0.0053 m. Considering the absolute value is with a relatively small magnitude, the difference that shifting in the direction along the slope length induced by the method is acceptable. In this sense, the method presented here can be applied to any object with similar characteristics with the rock slopes of narrow and long dimensions, e.g., metro tunnels, railways, highways, and large-scale bridges. The above-analysis indicates the main limitation of range image technique applied in the change analysis of rock slopes. However, it definitely has no influence on the cell change analysis, volume change estimation, etc. In our experiment, the width of the slope is 0.8 m, the range for Slope 1 is less than 9 m while for Slope 2 it is less than 15 m. The distance from the scanner center to the slope plane is less than 2 m. Here we define the distance OO equals to 2 m, and the range of 2 m to 15 m is used to illustrate the influence induced by the range image technique. A graph indicating the difference ( CD ) against the ranges is shown in Figure 18. It is apparently that the difference decreases as the range increase. The difference at the range of 2 m is 0.0396 m while the difference at the range of 15 m is 0.0053 m. Considering the absolute value is with a relatively small magnitude, the difference that shifting in the direction along the slope length induced by the method is acceptable. In this sense, the method presented here can be applied to any object with similar characteristics with the rock slopes of narrow and long dimensions, e.g., metro tunnels, railways, highways, and large-scale bridges. The above-analysis indicates the main limitation of range image technique applied in the change analysis of rock slopes. However, it definitely has no influence on the cell change analysis, volume change estimation, etc. sense, the method presented here can be applied to any object with similar characteristics with the rock slopes of narrow and long dimensions, e.g., metro tunnels, railways, highways, and large-scale bridges. The above-analysis indicates the main limitation of range image technique applied in the change analysis of rock slopes. However, it definitely has no influence on the cell change analysis, volume change estimation, etc. In order to investigate the influence induced by various point density on the volume change estimation results, the raw point clouds were downsampled with the percentage of 50%, 25%, and 10%, respectively. The surface density of slope 1 for the raw, 50%, 25%, and 10% of the point clouds is In order to investigate the influence induced by various point density on the volume change estimation results, the raw point clouds were downsampled with the percentage of 50%, 25%, and 10%, respectively. The surface density of slope 1 for the raw, 50%, 25%, and 10% of the point clouds is 378,089 pt/m 2 , 189,044 pt/m 2 , 94,522 pt/m 2 , and 37,809 pt/m 2 , respectively. While the corresponding surface density for slope 2 is 577,888 pt/m 2 , 288,944 pt/m 2 , 144,472 pt/m 2 , and 57,789 pt/m 2 , respectively. The proposed methodology was implemented on the downsampled point clouds again; the volume change information for the four cases are summarized in Table 6. The maximum volume changes for erosion occurred, extracted from the raw, 50%, 25%, and 10% of point clouds for slope 1 are 1257.1 mm 3 , 1261.9 mm 3 , 1192.0 mm 3 , and 1506.9 mm 3 , respectively. While for slope 2 the maximum volume changes are 1004.8 mm 3 , 1003.3 mm 3 , 994.4 mm 3 , and 991.7 mm 3 , respectively. The minimum volume changes, i.e., dilation occurred, extracted from the raw, 50%, 25%, and 10% of point clouds for slope 1 are −1229.8 mm 3 , −1194.2 mm 3 , −1262.5 mm 3 , and −1408.7 mm 3 , respectively, in which the symbol "−" represents the dilation occurred. While for slope 2 the minimum volume changes are −528.5 mm 3 , −540.7 mm 3 , −546.4 mm 3 , and −486.9 mm 3 , respectively. Afterwards, the volume changes for the entire damaged area were estimated with values of −79,060.9 mm 3 , −76,013.8 mm 3 , −86,922.1 mm 3 , and −81,596.1 mm 3 , for the raw, 50%, 25%, and 10% of the point clouds, respectively. Next, the comparison of the extracted volume change results between the downsampled point clouds and the raw point cloud was carried out. It is evident that the maximum and minimum volume changes have a relatively large difference for the downsampled point cloud of 10% w.r.t the raw point cloud. In the case of the downsampled point cloud of 10%, the average point space is~5.2 mm, namely on the surface of one stone (nominal diameter is 16.20 mm), the number of points is less than ten, so the points are not enough for reflecting the stone surface in an accurate way. Therefore, the low point density influences the extraction results of the volume change. In this regard, we conclude that relatively ideal volume change estimation results could be obtained under the surface point density greater than 144,472 pt/m 2 . To investigate the influence induced by various cell sizes on the volume change estimation results, the raster image is divided into cells of 0.01 rad × 0.01 rad, 0.005 rad × 0.005 rad, and 0.001 rad × 0.001 rad, respectively, for volume change estimation. The results are summarized in Table 7 (results of cell of 0.0025 rad × 0.0025 rad included). Thanks to the various cell sizes, it is reasonable that the maximum, minimum, and mean volume changes are different against the various cell sizes. The volume changes of the entire interest area, i.e., the sum of the volume changes w.r.t four cell sizes are expected to be identical. However, as shown in Figure 7, we found that for slope 1, the sum of the volume change becomes smaller with increased cell size, namely, the greater the cell size the smaller the sum of the volume change. While for slope 2, the sum of the volume change for the cell sizes of 0.01 rad × 0.01 rad, 0.005 rad × 0.005 rad, 0.0025 rad × 0.0025 rad, and 0.001 rad × 0.001 rad are −29,440.3 mm 3 , −23,144.4 mm 3 , −22,552.2 mm 3 , and −21,047.1 mm 3 , respectively, and their magnitudes are almost the same. As shown in Section 3.2, the mean distance from the damaged area to the scanner center for slope 1 is~6.2 m; while the value for slope 2 is~4.17 m. Therefore, for the raster image of slope 1 with the cell sizes of 0.001 rad × 0.001 rad, 0.0025 rad × 0.0025 rad, 0.005 rad × 0.005 rad, and 0.01 rad × 0.01 rad, the corresponding cell sizes in Euler space are~0.006 m, 0.015 m,~0.03 m, and~0.06 m, respectively. For the raster image of Slope 2 with the cell sizes of 0.001 rad × 0.001 rad, 0.0025 rad × 0.0025 rad, 0.005 rad × 0.005, rad and 0.01 rad × 0.01 rad, the corresponding cell sizes in Euler space are~0.004 m,~0.01 m,~0.02 m, and~0.04 m. At this point, we found that if the cell size in Euler space is greater than the nominal of the stone diameter, i.e., 16.2 mm, one cell covers more than one complete stone which causes the accretion or the erosion of the current cell averaged by the region outside the stone. In this regard, the effectiveness and reliable of the proposed methodology is constrained by the cell size, that is, it is better to set the cell size smaller than the nominal size of the particles composing the slopes in the real applications. As stated in a past paper [64], terrestrial laser scanning measurements are subject to noise induced from the incident angle. They concluded that above a 60 degree incidence angle, noise in the measurements affect the scan point precision. In our experiment, the incidence angle is close to 60 degrees, even greater than 60 degrees in some areas. However, considering the distinctiveness of the rock slope with the characteristics of long and narrow, it is unavoidable to guarantee the incident angles w.r.t entire areas greater than 60 degrees. As a result, we only used the damaged area that is close to the scanner center, i.e., the incidence angle is relatively small w.r.t the entire slope, for change analysis. This setting somehow improved the reliability of the volume change extraction results. Moreover, the raster image was introduced in this article which averaged all points in each cell. This step also reduced the error induced by the incidence angle. In the real application, it is better to mount the scanner with a relatively higher position with respect to the object to ensure a smaller incident angle which definitely increases the feasible and accuracy of the proposed methodology.
Also, the roughness of the slope is another factor that influences the plausibility of the proposed methodology. In our experiment, all the stones composing the slope are clear and have almost identical size, i.e., they have relatively low roughness. However, in the real application, the particles composing the slopes are dirty, e.g., with impurities, vegetation that grew around it and lichen, irregular in different areas. This maybe result in the roughness in different areas varies greatly, and then somehow influence the accuracy of the volume change estimation results. Therefore, an additional step is to assess the roughness of the object surface before implementing the proposed methodology. Several roughness assessment methods have been proposed, e.g., in multiple past papers [65,66], based on of which the raw point cloud could be cleared.

Discussion
The results obtained in Sections 3.1-3.5 show that the proposed methodology presented here is efficient and feasible. In this section, we emphasize the advantages of our approach, but also the practical challenges faced.
One of the main advantages of the method is the reduction of computational costs in an accurate way. In the conventional approach, as reported in a past paper [55], several additional processes, e.g., normal vector estimation, were required that had to executed on the full point cloud, before a suitable coordinate transformation could be estimated. The accuracy of these processing steps was affected by data noise, edge effects, variation in surface density, etc. In contrast, the range image technique, which is a 2.5D approach, enables the user to lift geometrical information back from 2D raster space to 3D space. This workflow, from 3D point cloud to 2D raster and back to 3D point cloud, ensures the directions of interest are easy and can be reliably estimated in 2D space. Besides, in a traditional Cartesian raster image based on coordinates transformation, cell sizes in 3D space are almost the same but the number of points in each cell varies with the local point density which depends on the measurement geometry. In the raster image approach, as discussed in this paper, the number of points in each 2D cell is almost identical which guarantees a certain uniformity of point cloud quality and accuracy.
The workflow in this research considered a specific indoor scene. Still, a similar workflow can and will also be applied in large outdoor coastal scenarios. For example, in [67], a so-called permanent laser scanning setup is discussed, in which the same sandy beach is scanned over and over again at hourly intervals for a total period of, for example, 6 months. To detect changes in such a spatiotemporal data set, similar methodology as presented in this paper will be applied. As scanning always takes place from the same location, only a fine registration procedure is needed to align scans from different epochs. Next, data of each epoch will be organized in an aspherical coordinate system as discussed in this research. Such organization will enable to handle the larger point clouds efficiently as neighbor detection can exploit the 2D grid organization, while the grid size can be adapted to size of the changes to be detected.
This methodology is notably suitable if scan data from one scan position is compared. Even if more scan positions are considered, a similar approach can be used, as already demonstrated by Zeibak and Filin in [68], where a second viewpoint is transformed to match a scan from a first viewpoint in a methodology that distinguishes areas that were changes from areas that are either stable or occluded because of the difference in viewpoints.
The change analysis results clearly indicate that the damage happened in the area near the intersection line between the wave surface and the slope plane. By predefining the size of the cells, we generated the raster image based on the range image technique to easily quantify the changes. Consequently, we performed the quantification and estimation of changes, which showed favorable results for detecting changes from rock slopes subject to wave attack test. To this end, the volume change estimation is done to validate the results. The results show that the transport direction actually swaps from seaward, for steeper 1:5 slope, to landward, for the less steep 1:10 slope. To this end, the limitation of the range image technique applied in change analysis of rock slopes is reported that related to the shifting in the direction along slope. Even so, after performing the range image technique, the computational performance is actually enhanced.
The methodology presented herein is of great practical importance for engineers involved in coastal zone management, coastal zone policy makers, etc., providing them guidance for the application of TLS in the evaluation of rock slopes. Specifically, the results about the transport direction of rock slopes with different steepness can serve to learn to understand the physical processes and can be further translated into governing practical engineering parameters, such as the selection of steepness for designing the coastal structures, damage evaluation, etc. Also, the possibility and effectiveness of the stones used in our experiment are verified. Furthermore, it has been proven that the laser scanning technique for understanding the motion of rock slopes can be extended to the change detection of coastal protection structures composing of stones.
However, for real applications, some problems need to be solved. First, in real scenarios such as shore lines, the region of interest is not necessarily narrow; the initial point cloud should be first subdivided into smaller clouds along the shore line with a defined width so that each partial cloud can be inspected for changes. A second problem could be the lack of identifiable points for TCP matching between two scans. Some stable artificial objects, e.g., mobile signal towers or corner of the buildings, could be used for registration between two scans.

Conclusions
The methodology proposed herein is an efficient and precise means of detecting changes from rock slopes using dense terrestrial laser scanning point clouds, providing a new way to evaluate the erosion of rock slopes to wave attack in coastal engineering. The estimation of erosion or dilation is a new possibility with the range image technique due to its large spatial resolution. In our case studies, the results are encouraging. The use of dense terrestrial laser scanning point cloud can provide information about the stone motion in the rock slopes more efficiently and at a higher resolution compared to classical instruments.
The processing of dense terrestrial laser scanning point cloud for change analysis presented in this article includes registration using TCP method, change detection using the cloud-to-cloud distances, range image generating based on the horizontal and vertical angular resolution, raster image creation, quantification of changes, and volume change estimation. During scanning, three plane targets were set between two scans before and after the wave attack test, a high precision registration was performed with the registration error less than 1 mm. The cloud-to-cloud distances were first adopted to realize the direct point cloud comparison, so that the damage area with significant change value was recognized. However, the notation of this type of distance is always positive that means it only purely reflects the magnitude of the change and the direction of the change is beyond expression. In such a case, for obtaining the magnitude and direction of the change, the range image technique was introduced. Based on the horizontal and vertical angular resolution given by the scanner itself, the transformation from 3D Cartesian coordinates to 2D range image is determined. For easier for quantification of changes, the raster image based on the range image technique was generated by predefining a cell size. As a result, the quantification and estimation of changes were performed and favorable results were also obtained for detecting changes from rock slopes subject to wave attack test.
The main contribution of this article focuses on the introduced range image technique for the change analysis of rock slopes with the characteristics of long and narrow, the limitations of which have been summarized in Section 3.5. However, from the perspective of the technique itself, it is really a promising approach for change analysis the long and narrow objects which somehow avoid introducing the unalignment error. Meanwhile, the identical number of points in different cells enables the equal-accuracy analysis which is a vital factor for the global analysis of the objects.
However, in our study, only two types of slopes were tested and analyzed, for insight into the influence induced by the various steepness, more types of steepness of slopes should be considered and implemented to obtain more information about the stability of slopes within various steepness. Meanwhile, more types of stones should be tested to find the optimal size of stone at a given steepness. In addition, due to the change that mostly happened near the intersection line, we believe that the change value is related the distance between the stone surface and the wave surface. However, we did not focus on this point in this paper which may be an interesting topic in the further study, the goal of which is to establish a rigorous mathematical model between the change value and the distance.
Upcoming research will focus on more detailed change information extraction such as each individual stone shape extraction and interpretation. However, it definitely a hard task from a laser scanning processing perspective due to the reasons such as the shape of the stone is irregular, it is difficult to acquire the point cloud for an entire stone, various point densities, etc. In reality, the movement of individual stones under wave attack is constrained in certain ways by gravity, wave forces, and friction between rocks; the motion can be regarded as a random process. For better understanding the change of the slope and the movement of individual stone, gravimetric analysis of each stone, friction analysis between stones, wave forces on the stones, and other possible influence factors are required. It is an interesting topic in the fields of surveying, remote sensing, coastal engineering, and hydraulic engineering.