Filtering Green Vegetation Out from Colored Point Clouds of Rocky Terrains Based on Various Vegetation Indices: Comparison of Simple Statistical Methods, Support Vector Machine, and Neural Network

Abstract

Filtering out vegetation from a point cloud based on color is only rarely used, largely due to the lack of knowledge of the suitability of input information (color, vegetation indices) and the thresholding methods. We have evaluated multiple vegetation indices (ExG, ExR, ExB, ExGr, GRVI, MGRVI, RGBVI, IKAW, VARI, CIVE, GLI, and VEG) and combined them with 10 methods of threshold determination based on training set selection (including machine learning methods) and the renowned Otsu’s method. All these combinations were applied to four clouds representing vegetated rocky terrain, and the results were compared. The ExG and GLI indices were generally the most suitable for this purpose, with the best F-scores of 97.7 and 95.4, respectively, and the best-balanced accuracies for the same combination of the method/vegetation index of 98.9 and 98.3%, respectively. Surprisingly, these best results were achieved using the simplest method of threshold determination, considering only a single class (vegetation) with a normal distribution. This algorithm outperformed all other methods, including those based on a support vector machine and a deep neural network. Thanks to its simplicity and ease of use (only several patches representing vegetation must be manually selected as a training set), this method can be recommended for vegetation removal from rocky and anthropogenic surfaces.

1. Introduction

Research on the complex morphology, genesis, and transformations of geological features cannot be conducted without accurate geodetic measurements that provide thorough surveys, maps, sections, and spatial models. Such objects can be studied using modern non-contact technologies such as digital photogrammetry [1] and laser scanning techniques [2], both ground-based [3] and airborne [4]. These methods usually produce a point cloud describing the measured surfaces in detail. Especially when mapping rock formations, landslides, etc., the presence of vegetation is one of the main obstacles [5,6,7] and, hence, must be filtered out before further processing. Various geometric and structural filters have been implemented in commercial or non-commercial software for use with LiDAR point clouds [8] or point clouds generated photogrammetrically [9]. New filtering algorithms or procedures have also been developed [10,11,12,13,14], some of which have been successfully applied to rock masses [15,16]. However, the successful use of algorithms with the geometric principle for rocky terrains is difficult due to the ruggedness of the formations.

To remove vegetation, color-based filtering (especially using green color) is a logical option. However, automated implementation of this solution in the RGB color space is complicated, although a human operator can distinguish between vegetation and ground under normal conditions without any problems. An alternative solution might lie in using vegetation indices.

A vegetation index is a numerical value that aims to express a characteristic of vegetation, such as its health. It can, therefore, be used to distinguish not only the quality of the vegetation but also whether it is vegetation or not. There are many vegetation indices calculated from different spectral bands of electromagnetic radiation. This approach is particularly used in conjunction with satellite imagery, where registered bands include invisible parts of electromagnetic radiation. However, for the application proposed here and in everyday practice, standard cameras are typically used due to the high costs of multispectral sensors and, in effect, only the visible bands captured by standard cameras (red, green, and blue) are usually available. This limits the range of usable indices, but a considerable number of such indices still remain available. Hereinafter, we will focus only on vegetation indices (VI) derived from the visible spectrum.

The use of vegetation indices to filter different vegetation types can be found in many studies. In a study by Meyer et al. [17], several vegetation indices were tested for the detection of soybean plants. Moorthy et al. in [18] detected sugar beet and corn using several established and one new vegetation indices. Kim et al. in [19] detected Chinese cabbage and white radish using multitemporal imagery and vegetation indices. Liu et al. in [20] worked with barley leaves and Ponti et al. in [21] used vegetation indices on images. All these studies show that the use of vegetation indices on image data is, to some extent, functional and useful. The use of color information in the cloud in conjunction with vegetation indices for vegetation filtering is addressed, e.g., by Anders et al. in [22]. Alba et al. in [23] very interestingly use first a near-infrared (NIR) camera cloud coloration from a laser scanner and only then apply an NDVI to filter the vegetation out. Cloud classification using vegetation indices for vine detection is addressed by Mesas-Carrascosa et al. in [24]. Núñez-Andrés in [25] compared vegetation filtering based on color information transformed into HSV (Hue, Saturation, and Value) and on ExG vegetation index on a rock massif, with both methods yielding good results. The use of colors for cloud processing can, therefore, greatly increase the efficiency of filtering or classification.

It is worth noting that although mass data collection methods are currently very popular in many fields of research involving complex terrain [26], such as coastal monitoring [27,28,29], volcano exploration [30], underground research [31,32], tree detection [33,34], the evolution of rock glaciers [35], monitoring of rock masses [36,37,38], or their geological analysis [39], vegetation filtering in such cases is typically handled by human operators as the performance of automated algorithms proposed for this purpose, so far, is generally poor as confirmed by Blanco et al. [6]. Specifically, they usually have problems when encountering highly rugged and/or sloped terrain [15]. Ignoring the vegetation and failing to remove it can greatly complicate the interpretation of the results [7], which further underlines the importance of introducing such an automated method.

One of the reasons why automatic vegetation filtering is problematic lies in the difficult determination of a threshold value for distinguishing green vegetation points from non-vegetation ones. Establishing any absolute and universally valid threshold is generally impossible because camera sensors are not identical and, in effect, such a threshold must necessarily differ among imagery acquisition systems. All the studies cited above use either Otsu’s method [40] or set a specific value for the threshold by estimation. Otsu’s method is an automated one, applicable to data where the representation of vegetation and ground classes are similar. It is, however, unsuitable for applications of general filtering of green vegetation from point clouds for the purpose of examining morphological formations, as the relative abundance of green vegetation in the cloud is highly variable (and in some cases very small), as stated in [41] or [42].

In view of the above, this paper aims to test different methods of threshold determination in conjunction with various visible vegetation indices and to find the most suitable method for filtering out green vegetation based on vegetation indices for rocky terrain. Our approach requires manual selection of sample areas (training data), which can overcome the problem with a variable amount of vegetation in the data. Testing will be performed on four real-world datasets capturing various types of rocks from different locations, providing a wide color range of rocks, vegetation, and surroundings. We employ multiple algorithms utilizing vegetation indices, from simple mathematical approaches, through the use of a support vector machine (SVM) approach, and up to a three-layer neural network (DNN). The latter two methods are used for individual vegetation indices as well as for their combinations. In addition, the results of these methods are compared with those acquired using the method devised by Otsu.

2. Materials and Methods

2.1. The Test Data

To evaluate whether green vegetation can be distinguished by the suggested approach, data on vegetated rock walls with varying color ranges were employed. Although the origin of the data plays no role in the evaluation, as we are looking for a universally applicable method of vegetation filtering, a brief description of the data collection process is provided below.

2.1.1. Data 1

The region of interest was in the town of Ledeč nad Sázavou (Central Bohemia, Czech Republic). The rock wall was approximately 55 m wide and 22 m high. The point cloud was acquired using the SfM-MVS method and the DJI Phantom 4 UAV camera. The flight was piloted manually, and imagery was acquired from a distance of 20 to 30 m from the terrain in a way that ascertains approximately 80% overlap. The original mean ground sampling distance (GSD) was 8 mm. For the purposes of this study, the cloud was diluted to an average point density of 1850 points/m² (i.e., the points are roughly in a 0.02 m grid). In total, the diluted cloud contained 2,861,035 points. The representation of green vegetation in the dataset identified manually by the operator was 19.5%. The point cloud is shown in Figure 1a.

2.1.2. Data 2

The point cloud was acquired using the SfM-MVS method with the DJI Phantom 4 UAV’s internal camera in the area of Horní Počenice near the capital city of Prague (Czech Republic). The selected rock wall is approximately 15 m wide and 4 m high. The flight was piloted manually, and imagery was acquired from a distance of 20 to 30 m from the terrain in a way that ascertains approximately 80% overlap. The original mean ground sampling distance (GSD) was about 7 mm. For the purposes of this study, the cloud was diluted to an average point density of 11,900 points/m² (i.e., the points are roughly in a 0.01 m grid). In total, the diluted cloud contained 1,973,242 points. The representation of green vegetation in the dataset identified manually by the operator was 38.5%. The point cloud is shown in Figure 1b.

2.1.3. Data 3

The point cloud was acquired using the SfM-MVS method and a DJI P1 camera carried by a DJI Matrice 300 UAV in the vicinity of the settlement of Dolní Kounice (approximately 10 km south-west of Brno, South Moravia, Czech Republic). The selected cliff is approximately 40 m in width and 16 m in height. The flight was piloted manually, and imagery was acquired from a distance of 20 to 30 m from the terrain in a way that ascertains approximately 80% overlap. The original mean ground sampling distance (GSD) was about 2 mm. For the purposes of this study, the cloud was diluted to an average point density of 1900 points/m² (i.e., the points are roughly in a 0.02 m grid). In total, the diluted cloud contained 1,425,119 points. The representation of green vegetation in the dataset identified manually by the operator was 11.4%. The point cloud is shown in Figure 1c.

2.1.4. Data 4

The point cloud was acquired using the SfM-MVS method and the DJI Phantom 4 UAV’s internal camera in Porta Bohemica (on the banks of the Elbe River, between Malé Žernoseky and Litochovice nad Labem, Northern Bohemia, Czech Republic). The selected area is approximately 40 m in width and 30 m in height. The flight was piloted manually, and imagery was acquired from a distance of 20 to 30 m from the terrain in a way that ascertains approximately 80% overlap. The original mean ground sampling distance (GSD) was approximately 7 mm. For the purposes of this study, the cloud was diluted to an average point density of 900 points/m² (i.e., the points are roughly in a 0.03 m grid). In total, the diluted cloud contained 1,504,239 points. The representation of green vegetation in the dataset identified manually by the operator was 29.6%. The point cloud is shown in Figure 1d.

2.2. Vegetation Indices Tested

The selected vegetation indices used for testing are listed in

Table 1. Vegetation indices used for testing.

Abbrev.	Name	Formulae	Reference
ExG	Excess Green	2g – r − b	[43]
ExR	Excess Red	(1.4R − G)/(R + G + B)	[17]
ExB	Excess Blue	(1.4B − G)/(R + G + B)	[44]
ExGr	Excess Green-Excess Red difference	E × G – E × R	[45]
GRVI	Green Red Vegetation Index	(G − R)/(G + R)	[46]
MGRVI	Modified Green Red Vegetation Index	(G2 − R2)/(G2 + R2)	[47]
RGBVI	Red Green Blue Vegetation Index	(G × G – R × B)/(G × G + B × R)	[47]
IKAW	Kawashima Index	(R − B)/(R + B)	[48]
VARI	Visible Atmospherically Resistant Index	(g − r)/(g + r − b)	[49]
CIVE	Color Index of Vegetation Extraction	0.441R − 0.811G + 0.385B + 18.787	[50]
GLI	Green Leaf Index	(2 × G – R − B)/(R + 2 × G + B)	[51]
VEG	Vegetative Index	g/((r0.667) × b0.333)	[52]

, along with the abbreviation, calculation formula, and literature reference detailing the original purpose for their creation or use. R, B, and G represent the digital number of red, green, and blue channels, respectively; the variables r, g, and b are defined as r = R/(R + G + B), g = G/(R + G + B), b = B/(R + G + B).

2.3. Methods of the Threshold Determination

2.3.1. Data Distribution and Training Datasets

When designing the method, it is necessary to assume that camera sensors, vegetation, lighting, and other measurement circumstances vary from case to case, and, therefore, it is necessary to determine the threshold that distinguishes the points representing green vegetation from the other points individually for each point cloud.

To show that for typical real-world data, such as ours, it is practically impossible to perform automatic filtering, and that it is necessary to use operator-defined subsets, we show a histogram for the whole Data 1 cloud for ExG vegetation index (a typical representative) in Figure 2g (histograms for all indices and Data 1 are given in Appendix B).

The histogram is represented by a smooth curve with no obvious peaks indicating (at least) two different groups of data with different VI magnitudes.

The ExG histogram has a heavier tail on the right side, i.e., the statistical probability distribution is clearly not normal, nor is the neighborhood of the main peak symmetrical. This is due to the unequal representation of colors in the data, but this is to be expected.

It was, therefore, necessary to prepare training classes containing all color shades typical of individual classes as demonstrated in Figure 2, in which the Data 1 dataset is used to illustrate the selection of green vegetation intended to be filtered out (Figure 2a), rock (Figure 2b), and three different colors of clay surfaces (Figure 2c–e). The histograms of these subclasses are shown in Figure 2h, indicating that the values of vegetation indices differ even between groups of terrain points and highlighting the importance of the selection of multiple terrain colors for correct threshold determination. It should be noted that in this process, selecting terrain points representing various colors is more difficult than selecting vegetation. In our case, the training dataset for each data consisted of approximately 10,000 points per dataset.

Figure 2h also illustrates that the data are continuous, i.e., we cannot see any obvious gaps between the VI values for vegetation and non-vegetation. The distributions, however, differ, which makes the use of standard statistical methods feasible for distinguishing between them.

Table 2. The overview of the methods for threshold determination evaluated in this paper, see Section 2.3 for a description of individual methods.

Abbreviation	Method Description
SCND	Single-class method based on the normal distribution assumption
SCHC	Single-class method based on histogram calculation
TCNDp	Two-class method based on the normal distribution assumption with a threshold separating the same quantile of both training classes
TCNDi	Two-class method based on the normal distribution assumption with a threshold in the intersection of normal distribution functions
TCHCp	Two-class method based on histogram calculation with threshold separating the same quantile of both training classes
TCHCi	Two-class method based on histogram calculation with a threshold in the intersection of smoothed histograms
TCSFf	Two-class method with a threshold maximizing the f-score function
TCSFs	Two-class method with a threshold determined based on the s-score function
SVM	Classification using the support vector machine (SVM)
DNN	Classification using the deep neural network
Otsu	Classification by the Ostu’s method applied on the whole point cloud

shows the overview of methods evaluated in this paper.

2.3.2. Single Class Method Based on Normal Distribution Assumption (SCND)

The single-class method assuming normal distribution (SCND) was the simplest method of threshold determination used in this paper. This method determines the threshold only from data representing the green vegetation class in the preselected training set (i.e., it ignores terrain data). The threshold determination is based on a simple assumption that 2.5% of the data at the margin of the distribution describing vegetation is already contaminated with terrain points. The threshold is then derived simply by calculating the mean and standard deviation (SD) of the distribution describing the vegetation and cutting off anything that lies further than 1.96·SD from the mean on the side adjacent to the histogram of terrain points (i.e., Figure 3a). The algorithm, therefore, works as follows:

• The mean (M) and standard deviation (SD) are determined from the VI values of the training set.
• The threshold T is calculated using the formula T = M + 1.96·SD or T = M − 1.96·SD; whether plus or minus is used depends on the orientation of the particular VI (if vegetation has higher values of VI than other points, the minus sign is applied, and vice versa).
• All points that exceed this threshold are removed from the cloud.

By observation, it was found that the histograms showing the distribution of the respective VIs for the vegetation class are very close to a normal distribution, and, therefore, the assumption of normality was met for our data. Another advantage is that, given the number of points in the point cloud, it is not a problem to select thousands to tens of thousands of points, and, thus, the statistical nature of the threshold is robust to small amounts of unwanted contamination (e.g., a few brown or black points remaining among the green points of a shrub).

2.3.3. Single Class Method Based on Histogram Calculation (SCHC)

The second simplest method, the single class method based on a histogram calculation (SCHC), is similar but does not assume a normal distribution (it is, therefore, universally applicable to any distribution). The threshold value is then determined as the value corresponding (in agreement with the previous method) to the 2.5-percentile from the histogram of the training vegetation class (see Figure 3b).

2.3.4. Two-Class Method Based on the Normal Distribution Assumption (TCND)

Classification methods based on the competition of two (or more) classes are generally considered more reliable. For this reason, the two-class method based on a normal distribution assumption (TCND) and the two-class method based on a histogram calculation (TCHC) were also tested. These methods do not determine the threshold using a fixed probability value, but attempt to separate two competing distributions.

The TCND method assumes that both the vegetation indices of the green vegetation and the rest of the cloud are normally distributed. It, therefore, uses both training classes as follows: the mean and standard deviation for the classes of green vegetation (M_V, S_V) and remaining points (M_R, S_R), respectively, are determined from the VI values of the training set.

Within this approach, two methods were employed for the threshold determination. In the percentile-based TCND method (TCNDp), the threshold was set to the points by removing identical percentiles from both vegetation and non-vegetation classes in the training set (see Figure 3c). The intersection-based methods (TCNDi) determine the threshold as the intersection of the two distributions (see Figure 3e).

The first method (TCNDp), therefore, calculates the threshold T using the formula:

$$T=\frac{M_V·S_R+M_R·S_V}{S_V+S_R}$$

(1)

This equation is based on the fact that at the threshold, the percentile cut off by the threshold is identical for both distributions and can be calculated as p-times the standard deviation (in SCND, p was set to 1.96). Therefore, the above equation is derived by equaling p_v = p_r from the equations p_v = (Mv − T)/Sv for vegetation and p_r = (T − M_R)/S_R for the remaining points; note that the signs in the brackets are inverse, since for one class, the cut-off percentile is to the right of the mean and for the other class, it is to the left. Nevertheless, whether the sign will be positive or negative in the individual equation depends on the vegetation index; see the explanation above in Section 2.3.1.

The second method (TCNDi) of two-class threshold determination is based on determining the intersection of the probability distribution functions of the two classes (i.e., green vegetation and remaining points). The calculation is based on solving the quadratic equation obtained by equaling two normal distribution functions:

$$f\left(x\right)=\frac{1}{S·\surd 2\pi }·e^{-\frac{{\left(x-M\right)}^2}{2S^2}}$$

(2)

$$\frac{1}{S_V·\sqrt{2\pi }}·e^{-\frac{{\left(x-M_V\right)}^2}{2S_V^2}}=\frac{1}{S_R·\sqrt{2\pi }}·e^{-\frac{{\left(x-M_R\right)}^2}{2S_R^2}}$$

(3)

The problem can be solved by logarithmization of the entire equation. The quadratic equation has generally two solutions—it is necessary to use the solution lying between M_V and M_R. It can be solved numerically as well.

2.3.5. Two-Class Method Based on Histogram Calculation (TCHC)

This method is basically a combination of the SCHC and TCND methods. The determination of the threshold is performed on the basis of shifting the provisional threshold until equal percentiles are cut off for both classes (TCHCp, see Figure 3e) or on the basis of the intersection of the histogram envelope curves (TCHCi, see Figure 3f).

We proposed numerical solutions by constructing a 1000-class histogram between the means of the vegetation (M_V) and remaining points (M_R), where the desired threshold must lie. The 1000 classes were used to create as fine a resolution as reasonably possible to improve the accuracy of the determined threshold. For the TCHCi method, it was necessary to smooth out the envelope curve of the histogram using a moving average as the representation within individual categories can vary; a window of 41 classes was used for this purpose.

2.3.6. Two-Class Methods Based on the Score Function Evaluation (TCSF)

The two-class method based on score function evaluation (TCSF) was another numerical method employed in this study. Two functions were used as the evaluation functions, namely f-score (TCSFf) and a custom function minimizing squares of the numbers of incorrect classifications (TCSFs).

TCSFf—the threshold is determined by maximizing the f-score. In practice, the score function is calculated for thresholds in regular steps (in our case, we used the step of a 1/10,000 of the interval between the means of the vegetation points and the remaining points); the point corresponding to the highest f-score was considered the threshold. The definition of the components for the equations below is given by the standard terminology of the binary classifiers (TP = true positive identification, FP = false positive, FN = false negative, TN = true negative) and is detailed in the next Section 2.4 since it is used for experimental results evaluation.

$${\mathrm{f}}_{\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}}=\frac{2\mathrm{T}\mathrm{P}}{2\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}+\mathrm{F}\mathrm{N}}$$

(4)

TCSFs = the threshold which is determined in a similar way based on the proposed s-score function describing the dependence of the percentage of unsuccessful classifications on the threshold (see Equation (1)). This approach supports a balanced minimization of the number of misclassified points in both classes (vegetation and remaining points).

$$\mathrm{s}=\frac{\sqrt{{\mathrm{F}\mathrm{P}}^2+{\mathrm{F}\mathrm{N}}^2}}{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}+\mathrm{F}\mathrm{N}}$$

(5)

2.3.7. Two-Class Method Based on the Support Vector Machine (SVM)

Automatic classification tools, such as support vector machines, can be used for distinguishing the training data classes without any additional assumptions. A support vector machine (SVM) is a supervised machine learning algorithm that can be used for both classification and regression. The goal of the SVM algorithm is to find a hyperplane in an N-dimensional space that distinctly classifies the data points. In our case, the freely available and easy-to-use LIBSVM toolbox (described in [53], downloadable at https://atoms.scilab.org/toolboxes/libsvm/1.5, accesses on 1 February 2023) for the Scilab universal computing system (version 6.1.1, www.scilab.org, accessed on 4 February 2023) was used. The input for training and classification can include virtually any number of parameters. The output of the classification is the value 0 or 1 (only two classes can be distinguished). In our case, only one input (vegetation index value) was used.

2.3.8. Two-Class Method Based on a Neural Network (DNN)

The use of machine learning in the form of a neural network is another option for data classification based on training sets. Classification using neural networks is a type of machine learning that involves training a model to recognize patterns in data and to make predictions about new data. Neural networks take inspiration from the learning process occurring in the human brain. Each element of the network (neuron) produces an output after receiving one or multiple inputs. Those outputs are then passed on to the next layer of neurons, which use them as inputs for their own function and produce further outputs. This continues until every layer of neurons has been considered and the terminal neurons have received their input. Those terminal neurons then output the final result for the model. The network used in our study was a three-layer network with one input neuron, 25 neurons of the hidden layer, and one output neuron. The number of neurons in the hidden layer may seem high; it was, however, determined through simple preliminary testing. We gradually increased their number to the point where the success rate of classifying the training data reached a plateau.

Similar to SVM, we used a toolbox for the Scilab system (ver. 6.1.1); this time, we employed the Neural Network Module (ver. 3.0, https://atoms.scilab.org/toolboxes/neuralnetwork/3.0/, accesses on 1 February 2023, created according to the book by Martin T. Hagan [54]). This toolbox supports several methods for learning the network. In our case, the Levenberg–Marquardt algorithm training function was clearly the best (based on preliminary testing, data not shown).

2.3.9. Two-Class Multi-VI Method Based on the Support Vector Machine (MSVM)

The SVM classification method can be easily used for a higher number of input parameters. We have, therefore, tested also the classification using a higher number of VIs, assuming a possible reduction in the uncertainty of the result. Based on the results of testing the previous methods, combinations of the 5, 3, and 2 most successful VIs were selected.

2.3.10. Two-Class Multi-VI Method Based on Neural Network (MDNN)

The reasoning is the same as in Section 2.3.8, as is the change in the number of input neurons to 5, 3, and 2, respectively.

2.3.11. Otsu’s Method (Otsu)

The method invented by Otsu [40] for automatic image thresholding was applied to the entire cloud as intended in their original paper. The algorithm returns a single threshold that separates values into two classes. The algorithm searches for the threshold that minimizes the intra-class variance, defined as a weighted sum of variances of the two classes. The class probability is computed from the bins of the histogram.

2.4. Testing of the Methods

For each test data, a reference dataset was manually created by classifying the whole point cloud into green vegetation and remaining (mostly terrain) data (all data shown in Appendix A). It should be noted that this is not only highly laborious but also partially subjective—especially in terms of assessing whether a particular point represents green vegetation or not. Here, the strategy used was that, if the assessment was uncertain, the point was not classified as a green vegetation point.

Due to the different percentages of points representing green vegetation in the whole cloud, it was necessary to manually prepare training classes containing an approximately equal representation of all color shades within the individual classes. These training sets were used to calculate the threshold for each VI and, based on the resulting thresholds, green vegetation was filtered out in the complete datasets and results were compared with the reference (manual) classification.

To determine the success of the filtering, established binary classification quality characteristics were calculated, namely the f-score (FS) and the balanced accuracy (BA which is clearly described, e.g., in the supplementary materials of [55]), the calculation is shown in

Table 3. Overview of success rate characteristics used.

Characteristics	Abbreviation	Calculation
F-score	FS	FS = 2TP/(2TP + FP + FN)
Balanced accuracy	BA	BA = (TPR + TNR)/2; TPR = TP/(TP + FN); TNR = TN/(TN + FP)

. Because this classification is binary, the data is classified into two categories: vegetation was designated as positive (P) and remaining points as negative (N). Successful classification is denoted as true (T), while an unsuccessful classification as false (F). Thus, TP denotes true positives, i.e., points correctly classified as green vegetation, FP denotes false positives, i.e., points incorrectly identified as green vegetation points, TN denotes points correctly classified other points, and so on.

3. Results

The above-described procedures were used to compare the results to those determined manually by a human operator. F-scores and balanced accuracies for each index and method of threshold determination are summarized in

Table 4. Results of the evaluation of the success of green vegetation filtering methods using individual VIs (F-score in %).

VI	SCND	SCHC	TCNDp	TCNDi	TCHCp	TCHCi	TCSFf	TCSFs	SVM	DNN	Otsu	Mean
ExG	97.7	92.6	90.8	91.7	93.5	95.1	94.4	94.7	94.0	94.1	85.5	93.1
ExR	69.1	74.5	72.5	72.1	73.7	76.0	75.3	75.2	73.2	74.9	64.6	72.8
ExB	68.6	72.7	79.2	81.2	77.0	78.0	77.9	77.8	78.8	78.2	81.3	77.3
ExGr	91.0	87.7	89.0	87.6	90.8	90.4	89.5	89.7	89.5	89.4	80.0	88.6
GRVI	80.6	80.9	77.3	76.5	78.7	81.1	80.1	79.7	77.9	79.9	68.9	78.3
MGRVI	81.1	81.0	77.5	76.8	78.7	80.8	80.1	79.7	78.6	80.0	71.0	78.7
RGBVI	90.1	92.2	91.1	91.8	91.1	92.4	91.1	91.3	91.0	91.2	86.8	90.9
IKAW	39.3	40.3	46.4	46.2	45.2	46.4	45.0	45.8	40.2	45.0	46.1	44.2
VARI	72.4	80.8	77.7	75.3	78.2	80.6	80.2	80.2	79.4	71.9	74.7	77.4
CIVE	87.9	87.4	90.4	90.3	89.2	89.9	89.6	89.6	89.9	89.5	84.8	89.0
GLI	95.4	92.6	92.0	92.4	93.5	94.9	94.4	94.7	94.0	94.2	86.0	93.1
VEG	70.4	91.6	85.3	87.3	93.9	92.7	92.7	92.9	92.8	92.8	77.2	88.2

and

Table 5. Results of the evaluation of the success of green vegetation filtering methods using individual VIs (Balanced accuracy in %).

VI	SCND	SCHC	TCNDp	TCNDi	TCHCp	TCHCi	TCSFf	TCSFs	SVM	DNN	Otsu	Mean
ExG	98.9	97.9	93.0	94.8	94.0	95.5	95.8	95.7	96.3	95.6	95.2	95.8
ExR	77.6	82.1	85.5	86.4	85.3	86.1	87.3	87.2	87.3	87.5	82.8	85.2
ExB	77.4	81.2	84.5	86.8	82.7	83.9	84.1	83.9	85.4	84.4	87.8	83.4
ExGr	96.3	95.8	94.2	95.0	94.9	95.8	95.7	95.8	96.2	95.8	91.7	95.6
GRVI	87.6	87.9	88.6	89.1	88.5	88.8	90.3	90.2	90.4	90.4	84.5	89.2
MGRVI	89.4	87.9	88.9	89.1	88.5	88.8	90.3	90.2	90.3	90.3	85.1	89.4
RGBVI	94.5	96.0	92.8	94.3	92.0	93.6	92.8	93.0	93.5	93.1	94.1	93.6
IKAW	49.2	49.9	61.0	60.1	60.4	60.8	59.7	60.0	48.4	50.4	61.0	56.0
VARI	83.6	87.6	88.5	88.8	88.4	88.8	90.3	90.2	90.2	81.0	87.1	87.7
CIVE	94.7	94.8	93.2	94.7	90.9	92.0	92.0	91.9	92.3	92.2	93.7	92.9
GLI	98.3	97.9	93.7	95.2	94.0	95.4	95.8	95.7	96.6	95.6	94.4	95.8
VEG	67.7	97.0	90.4	93.6	95.7	96.3	96.2	96.1	96.5	96.2	93.2	92.6

, respectively (values shown are mean results from all four datasets). On the far right in each table, the average value for the particular index across all methods is always shown. The highest (best) value in each column is highlighted in bold, showing the index producing the best results when using the respective method. Full details for each VI, data, and method are given in Appendix C.

In vegetation indices, the best average results were achieved when employing GLI and ExG. It is worth noting that the best results, namely an F-score of 97.7% and balanced accuracy of 98.9%, were acquired when using the simplest method, i.e., the SCND method, with the ExG vegetation index.

The ExG index generally performed best for SCND and SHC methods. When employing DND and DHC methods of filtering, the best results were achieved using GLI, followed by ExG.

From the machine learning methods, SVM and DNN yielded very similar results, with the best F-score and BA achieved in conjunction with GLI and ExG. Other VIs (in particular, IKAW, ExR, ExB, GRVI, MGRVI, and VARI with F-scores below 80%) were generally inferior.

Results of multi-VI methods are summarized in

Table 6. Results of the evaluation of green vegetation filtering success methods using multiple VIs.

	MSVM		MDNN		Mean
	FS [%]	BA [%]	FS [%]	BA [%]	FS [%]	BA [%]
ExG, ExGr, RBVI, CIVE, GLI	90.5	92.9	90.4	93.4	90.4	93.2
ExG, ExGr, RBVI, CIVE, VEG	90.6	93.0	89.6	92.8	90.1	92.9
ExG, RBVI, GLI	93.2	95.5	91.7	94.3	92.5	94.9
ExG, RBVI, CIVE	90.5	92.8	90.9	93.6	90.7	93.2
ExG, GLI, VEG	94.6	96.8	91.7	93.8	93.1	95.3
ExG, GLI	94.0	96.4	94.3	95.7	94.2	96.1

; in addition, the best results for several methods are illustrated in Figure 4. The most important observation here is that the results are very similar to the corresponding single VI method, suggesting that a well-chosen single VI, with a suitable method of threshold determination, which is significantly less computationally expensive, can be preferred to these relatively complicated methods.

4. Discussion

In this paper, we aimed to evaluate various vegetation indices and methods of threshold determination allowing the identification of green vegetation and filtering it out from a dense colored point cloud. The combinations of the indices and algorithms for threshold determinations were tested on datasets of rocky terrain containing 11.4–38.5% of green vegetation.

Methods that can be used to automatically find the threshold in a histogram were published previously. One such algorithm has been published by Otsu et al. [40] and has been widely used since (e.g., Refs. [17,18,19,21,24]). Kittler et al. concluded that Otsu’s method performs well when the histogram has a bimodal distribution with a deep and sharp valley between the two peaks [41]. Similar to all other global thresholding methods, Otsu’s method performs badly in cases of heavy noise, small object sizes, inhomogeneous lighting, and larger intra-class than inter-class variance [42]. The successful use of this algorithm, therefore, requires both classified classes (here green vegetation and other points) to be similarly represented—if the representation of green vegetation in the point cloud is small, they are practically invisible in the histogram and cannot be successfully identified. This is also the case with our data (see Figure 2g for ExG, or Appendix B, Figure A5, for all VIs—Data 1). In our data, the Otsu’s method indeed performed generally worse than the other tested algorithms, especially where the F-score is concerned. Figure 4b confirms the numerical findings, showing that even with the best-performing combination of Otsu’s method and vegetation index, the results are poorer than those acquired using the remaining methods. However, considering the above-mentioned limitations of the method, the result was not as poor as expected, especially where the balanced accuracy criterion is concerned.

Other authors set the threshold as an estimated fixed value. Such arbitrarily set threshold, however, differs among individual studies. For example, if considering only the ExG index, Anders et al. use a threshold value of 0.1 [22], while Núnez et al. used a threshold of −0.02 [25]. As each sensor is slightly different in terms of color calibration and as lighting conditions are also variable, it is not possible to determine a universally valid threshold; for example, the thresholds determined by the SCND method for our Data 1–4 in our study were 0.085; 0.167; 0.093; and 0.068, respectively. Similar variability was observed for other VIs as well.

Therefore, we have devised several semiautomatic methods that could be almost universally applied to determine the threshold and tested the success of classifications using individual methods. All tested methods were designed to overcome the problem of unbalanced representation of individual classes—manual selection of representative patches containing individual classes (training data) is the crucial step in our algorithm.

The results for individual indices and threshold determination methods shown in Table 4 and Table 5 imply several important findings of our study. First of all, all proposed methods of threshold determination provide good results when combined with appropriate vegetation indices. Even the algorithm with the poorest performance (TCNDp) yielded an F-score of 92% and BA of 94% when paired with the GLI. Interestingly, even the simplest algorithm, SCND, performed very well with all datasets. In combination with ExG, it yielded even the absolute best result in our study, classifying with an average F-score of 97.7% and BA of 98.9%.

On the other hand, the selection of the appropriate vegetation index plays a major role in the successful identification of vegetation. IKAW was shown to be unsuitable for this application, yielding the poorest results of all indices (F-scores of less than 50%). ExR, ExB, GRVI, MGRVI, and VARI indices also produced relatively poor results, with F-scores of approximately 80%, compared to over 90% for ExG, ExGr, RGBVI, CIVE, GLI, and VEG (depending on the method).

The good performance of SCND is an extremely valuable result, especially considering that this method is based on vegetation only and, therefore, does not require the operator to manually select non-vegetation training classes (i.e., surfaces of various hues/colors). This makes the entire process much simpler, faster, and more user-friendly and, together with low computational demands, makes this method, especially in combination with ExG or GLI, the option of choice. This method (Figure 4a) performed equally well or even better than much more complicated algorithms, including those based on neural networks (Figure 4d). Of course, the prerequisite of normal distribution of the vegetation index employed together with this method must be met to be able to use it; ExG and GLI, however, appear to generally meet this requirement well.

Where machine learning methods combining multiple indices are concerned, their results were also very good. However, as with individual indices, these methods did not bring any improvement over the simplest SCND method and, hence, the increased computational costs of their use appear (even though the increase in computational costs is not major as the training is performed only on a subset of the data) to be unjustified.

Comparing the performance of the F-score and balanced accuracy as the parameters for the evaluation of the algorithms, F-score appears to be more suitable for this particular application. When intending to filter out green vegetation, it is more important to remove all vegetation points even at the cost of removing a certain portion of terrain points than leaving vegetation in the cloud and considering it to be terrain. Looking at the equations for F-score and BA (Table 2), F-score puts more emphasis on the correct identification of the desired class while BA provides rather an estimation of the overall accuracy.

It is necessary to point out that the method proposed in this paper is not intended to replace standard vegetation filtering. On a grassy terrain, this method would remove virtually all points, which is, of course, not a required result. The strength of this method, however, lies in its excellent performance in identifying relatively scarce vegetation on monomaterial surfaces, such as (even highly rugged) rock formations or various anthropogenic structures. Further, filtering is not as perfect as if performed by a human operator and additional filtering by geometrical filters might be necessary (for example, to remove brown branches that remain unfiltered based on the color). Still, employing this filter first greatly simplifies the use of geometrical filters, which would likely fail on such a rugged terrain with vegetation cover [27]. When dealing with terrain suitable for the proposed algorithm, we, therefore, suggest using first the method proposed in this paper to remove green points (which, outside of areas of human habitation, certainly capture vegetation), and subsequently, employing commonly used geometric or structural filtering, or even manually removing the relatively few remaining vegetation points.

In addition, one must consider the acquisition period when intending to remove green vegetation—some vegetation changes color over the growing season and deciduous species drop their leaves, which may render the method inapplicable in certain seasons.

We have strived to minimize any possible limitations in our study by using multiple indices, multiple methods, and two evaluation criteria. However, one limitation is inherent to any similar study—the reference data must be manually classified, which may always bring about a certain amount of error, especially on the edges of the vegetated areas where the gradual change of green to the terrain color may be problematic. However, the number of points that might be misclassified in this way in our areas is relatively low as the vegetated areas are rather continuous patches than dissipated small individual plants.

5. Conclusions

In this paper, a color-based method for filtering green vegetation out from point clouds with scarce vegetation was proposed, multiple algorithms for determining a threshold for vegetation removal were evaluated, and the results were compared. All evaluated methods of threshold determination performed relatively well when used together with appropriate vegetation indices. In general, ExG and GLI indices appear to be the most suitable for this purpose. Surprisingly, the simplest method (single-class normal distribution, SCND) required only a selection of the sample of the vegetated area as the training set performed best, outperforming even support vector machine or deep neural network approaches. This is all the more valuable that this method is fast, is not demanding complicated calculations, and is easy to use, i.e., not necessitating the selection of terrain samples as well. This rapid and simple method is, therefore, highly suitable for filtering out vegetation, for example, from rocky or anthropogenic terrain.