-
-
Notifications
You must be signed in to change notification settings - Fork 10
/
gnd.qmd
144 lines (99 loc) · 8.21 KB
/
gnd.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
```{r setup,echo=FALSE,message=FALSE,warning=FALSE}
library(lidR)
library(ggplot2)
source("function_plot_crossection.R")
r3dDefaults = rgl::r3dDefaults
m = structure(c(0.921, -0.146, 0.362, 0, 0.386, 0.482, -0.787, 0,
-0.06, 0.864, 0.5, 0, 0, 0, 0, 1), .Dim = c(4L, 4L))
r3dDefaults$FOV = 50
r3dDefaults$userMatrix = m
r3dDefaults$zoom = 0.75
knitr::opts_chunk$set(
comment = "#>",
collapse = TRUE,
fig.align = "center")
rgl::setupKnitr(autoprint = TRUE)
knitr::opts_chunk$set(echo = TRUE)
```
# Ground Classification {#sec-gnd}
Classification of ground points is an important step in processing point cloud data. Distinguishing between ground and non-ground points allows creation of a continuous model of terrain elevation (see @sec-dtm)). Many algorithms have been reported in the literature and `lidR` currently provides three of them: Progressive Morphological Filter (PMF), Cloth Simulation Function (CSF) and Multiscale Curvature Classification (MCC) usable with the function `classify_ground()`.
## Progressive Morphological Filter {#sec-pmf}
The implementation of PMF algorithm in `lidR` is based on the method described in [Zhang et al. (2003) ](https://ieeexplore.ieee.org/document/1202973) with some technical modifications. The original method is raster-based, while `lidR` performs point-based morphological operations because `lidR` is a point cloud oriented software. The main step of the methods are summarised in the figure below:
<center>
</center>
The `pmf()` function requires defining the following input parameters: `ws` (window size or sequence of window sizes), and `th` (threshold size or sequence of threshold heights). More experienced users may experiment with these parameters to achieve best classification accuracy, however `lidR` contains `util_makeZhangParam()` function that includes the default parameter values described in [Zhang et al. (2003) ](https://ieeexplore.ieee.org/document/1202973).
```{r, warning=FALSE}
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzrn")
las <- classify_ground(las, algorithm = pmf(ws = 5, th = 3))
```
We can now visualize the result:
```{r plot-gnd-classification, rgl=TRUE}
plot(las, color = "Classification", size = 3, bg = "white")
```
To better illustrate the classification results we can generate and plot a cross section of the point cloud (see @sec-plot-crossection)).
```{r plot-gnd-crossection, echo=FALSE, fig.height=2, fig.width=8}
p1 <- c(273420, 5274455)
p2 <- c(273570, 5274460)
plot_crossection(las, p1 , p2, colour_by = factor(Classification))
```
We can see that although the classification worked, there are multiple points above terrain that are classified `2` (i.e. "ground" according to [ASPRS specifications](http://www.asprs.org/wp-content/uploads/2019/07/LAS_1_4_r15.pdf)). This clearly indicates that additional filtering steps are needed and that both `ws` and `th` parameters should be adjusted. Below we use multiple values for the two parameters instead of a single value in the example above.
```{r, message = FALSE}
ws <- seq(3, 12, 3)
th <- seq(0.1, 1.5, length.out = length(ws))
las <- classify_ground(las, algorithm = pmf(ws = ws, th = th))
```
After this adjustment the classification result changed, and points in the canopy are no longer classified as "ground".
```{r plot-gnd-crossection-2, echo=FALSE, fig.height=2, fig.width=8}
plot_crossection(las, p1 = p1, p2 = p2, colour_by = factor(Classification))
```
## Cloth Simulation Function {#sec-csf}
Cloth simulation filtering (CSF) uses the [Zhang et al 2016](http://www.mdpi.com/2072-4292/8/6/501/htm) algorithm and consists of simulating a piece of cloth draped over a reversed point cloud. In this method the point cloud is turned upside down and then a cloth is dropped on the inverted surface. Ground points are determined by analyzing the interactions between the nodes of the cloth and the inverted surface. The cloth simulation itself is based on a grid that consists of particles with mass and interconnections that together determine the three-dimensional position and shape of the cloth.
<center>
</center>
The `csf()` functions use the default values proposed by [Zhang et al 2016](http://www.mdpi.com/2072-4292/8/6/501/htm) and can be used without providing any arguments.
```{r, message = FALSE}
las <- classify_ground(las, algorithm = csf())
```
Similar to the previous examples, classification results can be assessed using a cross section:
```{r plot-gnd-crossection-3, echo=FALSE, fig.height=2, fig.width=8}
plot_crossection(las, p1 = p1, p2 = p2, colour_by = factor(Classification))
```
While the default parameters of `csf()` are designed to be universal and provide accurate classification results, according to the original paper, it's apparent that the algorithm did not work properly in our example because a significant portion of points located in the ground were not classified. In such cases the algorithm parameters need to be tuned to improve the result. For this particular data set a set of parameters that resulted in an improved classification result were formulated as follows:
```{r plot-gnd-crossection-4, message=FALSE}
mycsf <- csf(sloop_smooth = TRUE, class_threshold = 1, cloth_resolution = 1, time_step = 1)
las <- classify_ground(las, mycsf)
```
```{r, echo=FALSE, fig.height=2, fig.width=8}
plot_crossection(las, p1 = p1, p2 = p2, colour_by = factor(Classification))
```
We can also subset only the ground points to display the results in 3D
```{r plot-gnd, rgl = TRUE}
gnd <- filter_ground(las)
plot(gnd, size = 3, bg = "white")
```
## Multiscale Curvature Classification (MCC) {#sec-mcc}
Multiscale Curvature Classification (MCC) uses the [Evans and Hudak 2016](https://ieeexplore.ieee.org/document/4137852) algorithm originally implemented in the [mcc-lidar](https://sourceforge.net/projects/mcclidar/) software.
```{r, echo=FALSE}
las$Classification <- 0L
```
```{r plot-gnd-crossection-5}
las <- classify_ground(las, mcc(1.5,0.3))
```
```{r, echo=FALSE, fig.height=2, fig.width=8}
plot_crossection(las, p1 = p1, p2 = p2, colour_by = factor(Classification))
```
## Edge Artifacts {#sec-edge-artifact}
No matter which algorithm is used in `lidR` or other software, ground classification will be weaker at the edges of point clouds as the algorithm must analyze the local neighbourhood (which is missing on edges). To find ground points, an algorithm need to analyze the local neighborhood or local context that is missing at edge areas. When processing point clouds it's important to always consider a buffer around the region of interest to avoid edge artifacts. `lidR` has tools to manage buffered tiles and this advanced use of the package will be covered in @sec-engine).
## How to Choose a Method and Its Parameters? {#sec-method-selection}
Identifying the optimal algorithm parameters is not a trivial task and often requires several trial runs. `lidR` proposes several algorithms and may introduce even more in future versions; however, the main goal is to provide a means to compare outputs. We don’t know which one is better or which parameters best suit a given terrain. It’s likely that parameters need to be dynamically adjusted to the local context, as parameters that work well in one file may yield inadequate results in another.
:::{.callout-note}
`lidR` is an research and development package, not a production package. The goal of lidR is to provide numerous and handy ways to manipulate, test, compare point cloud processing methods. If available, we recommend using classifications provided by the data provider. The `classify_ground()` function is useful for small to medium-sized unclassified regions of interest because it is feasible to visually assess classification results. For large acquisitions where visual assessment is no longer feasible, we do not recommend performing ground classification without first studying its accuracy.
:::
## Other Methods {#sec-other}
Ground segmentation is not limited to the 3 methods described above. Many more have been presented and described in the literature. In @sec-plugins we will learn how to create a plugin algorithm and test it seamlessly in R with `lidR` using a syntax like:
```r
las <- classify_ground(las, algorithm = my_new_method(param1 = 0.05))
```