Comparison of Filters for Archaeology-Specific Ground Extraction from Airborne LiDAR Point Clouds

Abstract

Identifying bare-earth or ground returns within point cloud data is a crucially important process for archaeologists who use airborne LiDAR data, yet there has thus far been very little comparative assessment of the available archaeology-specific methods and their usefulness for archaeological applications. This article aims to provide an archaeology-specific comparison of filters for ground extraction from airborne LiDAR point clouds. The qualitative and quantitative comparison of the data from four archaeological sites from Austria, Slovenia, and Spain should also be relevant to other disciplines that use visualized airborne LiDAR data. We have compared nine filters implemented in free or low-cost off-the-shelf software, six of which are evaluated in this way for the first time. The results of the qualitative and quantitative comparison are not directly analogous, and no filter is outstanding compared to the others. However, the results are directly transferable to real-world problem-solving: Which filter works best for a given combination of data density, landscape type, and type of archaeological features? In general, progressive TIN (software: lasground_new) and a hybrid (software: Global Mapper) commercial filter are consistently among the best, followed by an open source slope-based one (software: Whitebox GAT). The ability of the free multiscale curvature classification filter (software: MCC-LIDAR) to remove vegetation is also commendable. Notably, our findings show that filters based on an older generation of algorithms consistently outperform newer filtering techniques. This is a reminder of the indirect path from publishing an algorithm to filter implementation in software.

1. Introduction

Airborne light detection and ranging (LiDAR) is now widely used in archaeological prospection [1,2]. Archaeologists interpret enhanced visualizations of airborne LiDAR-derived high-resolution digital terrain models (DTMs) with a combination of perception and comprehension, e.g., [3,4]. The results have proven to be very successful in detecting archaeological features and have uncovered numerable new settlements and almost innumerable new archaeological features worldwide, e.g., recently [5,6,7,8,9].

Airborne LiDAR-derived data provide archaeologists with detailed topography, which is a powerful resource for interpretation [10]. The often described process from data acquisition to archaeological interpretation [11,12,13,14,15,16] can be divided for the purposes of this article into raw data processing (project planning, system calibration, data acquisition, georeferencing, flight strip adjustment), point cloud processing (ground point filtering, additional processing), creation of end products (archaeology-specific raster elevation model interpolation and it visualization), and archaeological interpretation (interpretative mapping, ground inspection, deep integrated multi-scale interpretation). For archaeology, the most critical part of the workflow is ground point filtering [16], and this is the subject of our article.

It is worth reiterating that this process involves assumptions and decisions during the entire workflow, but these are often underestimated or even ignored [12,13]. As some of the end products (for example, a simple hill-shaded representation of the DTM) are almost self-explanatory, this may lead to a certain ignorance of the underlying data processing strategies [16].

In comparison to other fields the use of airborne LiDAR-derived data in archaeology is specific in several ways. First, visual inspection of enhanced raster visualizations remains the key method. The digital surface model (DSM) used often combines a digital terrain model (DTM) with off-terrain archaeological features [16] and modern buildings [10,17,18]. Such a specific DSM has been termed the digital feature model (DFM) [19].

Secondly, morphologically speaking archaeological features are anomalies, i.e., sporadically occurring atypical features of the terrain. The direct application of existing generic methods developed to assess the average accuracy of the terrain [20] is, therefore, less than ideal. For example, in archaeology, low noise can be tolerated as long as the size of the distortion is significantly smaller than the archaeological features.

Third, the time, effort, equipment, and human resources invested in airborne LiDAR data processing are only a small part of a typical archaeological project. For example, a detailed desk-based assessment of 1 km² takes about 1 person/day, and even the most rudimentary field survey of the same area takes up to 5 person/days. It is therefore less important whether the filtering of the point cloud takes 4 or 400 s than whether a single feature is incorrectly represented. For this purpose, archaeology-specific workflows have been proposed for both custom [16] and general-purpose data [15], none of which are as yet generally accepted. Two reasons for this are that the proposed workflows are too complex and in parts based on costly software.

Finally, we are currently experiencing an unprecedented expansion of archaeological applications of airborne LiDAR data due to the recent widespread availability of low or medium density, nationwide, general-purpose data, especially in Europe and the USA. However, the acquisition and processing of this data has not been optimized for archaeology. In most cases, the data are suitable for (at least some) archaeological purposes if custom archaeology-specific processing is applied, e.g., [15,21].

From the above, it can be concluded that archaeology as a discipline will greatly benefit from a workflow for the archaeology-specific processing of low and medium density airborne LiDAR data based on the free or low-cost off-the-shelf software. To achieve this, an archaeology-specific comparison of filters for ground extraction from airborne LiDAR point cloud is a necessary first step.

For this purpose, we compare the performance of nine automatic filters implemented in free or low-cost off-the-shelf software. The goal of this comparison is three-fold: to determine the archaeology-specific performance of existing filters, to determine the influence of the point density on the filter performance, and to identify directions for future archaeology-specific research.

In the past several studies have compared relevant filters implemented in free software [22,23,24,25,26,27]. However, none of them include newer filters, none of them use a standardized comparison method by Sithole and Vosselman [20], and none of them is archaeology-specific. Especially the former makes our comparison timely and relevant beyond archaeology.

A note on terminology is appropriate at this point. In general, we have followed the terminology widely accepted in relevant scientific publications. We use the term algorithm to refer to theory (e.g., academic articles) and filter to implementation (e.g., in software). The term filtering is used to describe the process of applying the filter to the data. The ground extraction filters are used to classify each point in the point cloud as either ground or non-ground and no point is deleted (filtered). However, when the DTM is interpolated from the point cloud, only the points classified as the ground are used and the rest are omitted (filtered) from the process, cf. [20]. This leads to a certain ambiguity in the terminology, especially concerning the term filter/filtering.

2. Materials and Methods

2.1. Data

Test sites from Austria (AT), Slovenia (SI), and Spain (ES) were selected based on their archaeological and morphological similarity to compare the filter performance at different point densities. Each site features:

• Hilltop settlement (archaeology)
• Buildings (modern)
• Vegetation on steep slopes
• Sharp discontinuities.

The additional test site (SI2) was chosen for extremely dense vegetation as an example of data with significant gaps in the coverage of the ground points. Each site is 500 × 500 m, but only the most relevant windows are shown in illustrations.

AT airborne LiDAR data was acquired in 2009. These data have a nominal average density of 4 pnt/m² and an estimated horizontal and vertical root mean square error (RMSE) of 0.15 and 0.40 m, respectively. Ground point filtering is not documented [28]. Point cloud data are available via a formal request (we have acquired the data in the non-classified XYZ format, other formats may be available). The AT test site has 3,844,873 points with an average density of 18.79 pts/m², which is significantly higher than the nominal average density.

The SI data was acquired in 2014. These data have a nominal density of 5 pnt/m² and an estimated horizontal and vertical RMSE of 0.09 m. Ground point filtering was performed with the proprietary filter gLidar [29] and is distributed via the eVode WMS (LAS or LAZ format, ISPRS classes 0–7 are populated). The SI1 test site has 2,809,344 points with an average density of 12.28 pts/m² and the SI2 has 2,756,756 points with an average density of 12.43 pts/m². It must be mentioned that due to the specific processing of the raw data, each point has a near-double “shadow point”. The “shadow point” is within 0.14 m horizontally and within 0.03 m vertically from the original point, but the distance and bearing to the original is not constant. This means that the functional point density is half of the above-mentioned, namely 6.14 pts/m².

The ES data was acquired in 2015, with an average density of 1 pnt/m² and the altimetric accuracy and precision with a RMSE of less than 0.20 m. Ground point filtering is not documented [30]. The point cloud data is distributed by the CNIG Download Center (LAZ format, ISPRS classes 0–7, 10, 12 are populated). The test site ES has 453,255 points with an average density of 1.83 pts/m².

The most important metadata are summarized in

Table 1. Metadata on test sites.

	AT	SI1	SI2	ES
Pnt/m2	18.79	12.28	12.43	1.83
File format	XYZ	LAS/LAZ	LAS/LAZ	LAZ
ISPR classes	/	0–7	0–7	0–7, 10, 12
Year of acquisition	2009	2014	2014	2015

.

It is not our intention to present the archaeology of the test sites in any detail, but merely to provide the key reference for each archaeological site and to illustrate the archaeological features (Figure 1).

The test site AT is located south of Graz near the town Wildon on a rocky hill consisting of algal limestone and marlstone interspersed with layers of clay that rise above the Mur River valley [31]. It is covered with a dense, mature deciduous forest with little undergrowth. The multiperiod archaeological site Wildoner Schlossberg (N 46°53′05″, E 15°30′48″) is mainly known for two medieval castles and a settlement [32].

The neighboring test sites SI1 and SI2 are located in the Dinaric Karst region named after the Pivka River southwest of Ljubljana. Both sites are located on limestone formations from the Jurassic period, which rise above the karst field Pivška kotlina, known for its intermittent lakes [33]. SI1 is covered by decimated deciduous forest in the first stage of regrowth. It is the site of the prehistoric hillfort Kerin above Pivka (N 45°40′21″, 14°11′40″) [34]. SI2 is characterized by several meadows, each surrounded by very dense deciduous hedges and it is the site of the Roman period villa rustica Sela near Pivka (N 45°40′56″, E 14°12′25″).

Test site ES is located on the magmatic granite hillock midway between the valleys of the Ulla River and its tributary the Vea south of Santiago de Compostela in northwest Spain [35]. It is covered with open deciduous forest and with two small clearings to the east of the site. The archaeological site O Castelo (N 42°44′30″, E 8°33′02″) is a Roman military camp with a probable Iron Age predecessor [36].

On the above mentioned sites 3 out of 4 morphological types of archaeological features are distinguished:

• features embedded in the ground (slight positive or negative bulge, typically up to 0.5 m rise over 5 to 20 m run; blue in Figure 1),
• features partially embedded in the ground (positive or negative spike, typically more than 0.5 m rise over 5 m run; green in Figure 1),
• standing features (an archaeological term for an off-terrain object characterized by a sharp discontinuity in the ground; red in Figure 1), and
• large standing objects (large building-like structures).

Most of the archaeological features detected on airborne LiDAR data fall into the first two types. Examples of type 1 are remains of trenches, ditches, fossil fields, past land divisions, and ancient roads. Examples of type 2 features are dwellings, ramparts, terraces, and burial mounds. Type 3 features occur relatively rarely, mostly as isolated sections of walls (without roof). From the perspective of point cloud filtering, they are non-ground objects. Examples of type 4 features include Mayan monumental architecture at Aguada Fénix, Mexico [5] and Khmer temples in Angkor, Cambodia [7]. These are region specific and are, to our knowledge, not present in European archaeology. Type 4 features are not the subject of this article since they demand additional processing steps that are not accomplished with ground extraction filters.

2.2. Evaluated Filters

We compare all filters known to us that are implemented in an off-the-shelf free or low-cost software under active development or maintenance (

Table 2. Evaluated filters and software used in the study (for more details see Appendix A).

Filter	Acronym	Software	Software Type
progressive morphological f.	PMF	lidR 2.2.4	free
slope-based f.	SBF	Whitebox tools	free
simple morphological f.	SMRF	PDAL 2.1.x	free
progressive triangulated irregular network	PTIN	lasground_new	commercial
weighted linear least-squares interpolation	WLS	Fusion 3.80	free
cloth simulation f.	CSF	CloudCompare 2.10.2	free
segmentation based f.	SegBF	Whitebox tools	free
multiscale curvature classification	MCC	MCC-LIDAR 2.2	free
Blue Marble Geographic’s hybrid f.	BMHF	Global Mapper 21.x	commercial

; Appendix A). Accordingly, BCAL Lidar Tools [37] was not considered because its freeware standalone version has been deprecated. Also not considered were legacy Whitebox GAT [38], ALDPAT [39], and LASground [40]. gLiDAR [41,42] was rejected because it was not designed for production and is arguably also legacy program. Furthermore, the skewness balancing filter [43], as implemented in PDAL, was not included because it does not work as intended (in our attempts, all points below a certain height above sea level have been marked as ground and all other as non-ground).

Archaeologists do not typically work with full-waveform airborne LiDAR data (but see [44]), therefore we do not assess post-processing methods for this specific category of airborne LiDAR data in the context of this article.

The analyzed filters are summarized in Table 2.

To our knowledge PMF, SBF, SMRF, CSF, SegBF, and BMHF were not assessed in any of the previous studies [22,23,24,25,26,27].

These and numerous other existing filters are based on algorithms that can be divided into five categories [26,45,46,47]:

• morphological filtering (PMF, SBF, SMRF),
• progressive densification (PTIN),
• surface-based filtering (WLS, CSF),
• segmentation-based filtering (SegBF),
• other (MCC), and
• hybrid (BMHF).

Morphological filtering algorithms are based on the concept of mathematical morphology. They retain terrain details by comparing the elevation of neighboring points [46,47]. A point is classified as ground if its height does not exceed the height of the eroded surface. The filtering is done by two elementary operations: a morphological erosion with the kernel function and a test for the difference between the original height and the eroded height of a point [48].

Progressive densification algorithms work progressively and iteratively classify more and more points as ground by comparing the elevation of the points with the estimated values by different interpolation methods. The initial set of seed ground points is often obtained from the lowest points in a large grid. Ground points are then densified in an iterative process, in which within every triangle the offsets of the candidate points to the triangle surface are analyzed. Points with offsets below the threshold are classified as ground and the algorithm proceeds with the next triangle. The new points are inserted into the triangulation. The process is repeated on the thus densified network until (i) there are no more points in a triangle, (ii) a certain density of ground points is reached, or (iii) when all acceptable points are closer to the surface than a threshold value [26,46,47].

Surface-based filtering algorithms first assume that all points are ground and then gradually remove the points that are not ground. Typically, simple kriging is used in the first step to describe the surface using all points. An averaged surface is then created between the ground and non-ground points. The residual values are created according to their distance from the averaged surface and the points are weighted according to their residual values. Points that have negative residual values are weighted more heavily and are considered ground points. There are many extensions and variants of the surface-based filters [26,46].

Segmentation-based filtering algorithms differ from the above-mentioned ones in that they first form homogeneous segments and then analyze segments rather than individual points. Segmentation can be performed using region-growing techniques, adopting either the normal vector or its change, resulting in either planar or smoothly varying surfaces, respectively. Alternatively, homogeneity can be achieved from height or height change only. The segmentation process has been described by the rationale that every point of a segment can be reached from any other point of the same segment, in the sense that it would be possible for a person to walk along such a path. Many implementations rasterize the data in the process and the segmentation methods are often derived from image processing [46].

Other approaches have also been developed. Some are based on determining a set of features for each point depending on the distribution of neighboring points or as measured features of each point. A decision tree based on machine learning is then often used to classify points based on their feature values. Some methods are based on the slopes or curvatures between points [46].

Hybrid filters combine algorithms to merge the strengths of different approaches [46,47].

The filters used in this article and their specific implementation in the software are described in more detail in Appendix A.

When comparing the different categories of algorithms, the following points can be summarized. Morphological filtering invariably impose a tradeoff. The erosion of the ground on steep slopes on the one hand and the inclusion of the non-ground points in flat areas on the other hand has to be balanced on a case-by-case basis. A general advantage of both progressive densification and surface-based algorithms is that the geometric properties of the terrain surface and the point cloud data characteristics are separated to a certain extent. Surface-based algorithms produce more reliable results since they also consider the surface trend, and due to their iterative nature, allow for a better adaptation to the terrain. The segmentation-based algorithms have advantages mostly in urban areas where flat segments are prevailing [46]. Bujan recommends that morphological filtering algorithms are best suited for terrains with small objects, progressive densification for urban areas and terrains with diverse objects, while surface- and segmentation-based algorithms work best in forests and on steep slopes [47]. To this, we may add that one of the key features of segmentation-based algorithms is that the outliers (low and high) theoretically have a limited impact on the result. Thus, they should be suitable to handle data sets with significant noise such as point clouds obtained by structure from motion (SfM).

Nevertheless, all algorithms present similar difficulties when implemented in sharply changing terrain, low vegetation, areas with dense non-ground features, and complicated landscapes [20,46,47,49,50]. The two-decade-old expectation that the fusion of multiple sources (e.g., RGB values or intensity) and integration of different methods will improve the performance [20,51,52] is yet to materialize in off-the-shelf software.

The long-term comparison of the performance of specific algorithms is also reveling. Between 1996 and 2007, uniform algorithms prevailed, with densification being the most important one. Between 2007 and 2012, hybrid algorithms prevailed, and between 2012 and 2016 there was a period of renewed interest in uniform algorithms [47]. Early studies showed that algorithms incorporating a concept of surface performed better than morphological algorithms. In the test environment, a gradual improvement in the performance of algorithms, regardless of the concept they are based on, has been observed in the past two decades [47].

2.3. Method for Quantitative and Qualitative Assessment

Above mentioned studies that have compared relevant filters [22,23,24,25,26,27] used different approaches to calibrating filters with custom settings, which depends strongly on the skillset of the operator. Some omitted calibration by using default settings to minimize the subjective influence of the operator [25]. However, not only are the results more accurate when filters are calibrated, but some filters are also more susceptible to calibration than others [27]. Accordingly, most studies compare calibrated filters [22,23,26]. Our decision to compare calibrated filters was ultimately guided by our search for the best workflow, which must always include calibration.

We have minimized the subjectivity of the calibration process through intensive testing (383 preliminary DTMs were created as a result). For most filters, we used the method of relevant incremental adjustments for each variable, starting with the default value until the best result was found. As the process was guided by experience, in most cases, fewer than 4 attempts were required for each variable. As a typical example is the variable(s) of the MCC, case study AT:

• step 1: default value 2.0;
• step 2: value 3.0 leads to a deterioration;
• step 3: value 1.0 leads to an improvement;
• step 4: value 0.5 leads to a deterioration;
• result: value 1.0 is optimal.

Such an approach was not possible with the SMRF, since out of five variables there are two inversely related pairs. This means that a pair of variables must be changed each time, resulting in almost infinite combinations. Thus, we tested 15 settings that were optimized for ISPRS reference data sets [53].

The 36 best results—9 filters each for AT (Figure 2), SI1 and SI2 (Figure 3), and ES (Figure 4)—were compared qualitatively and quantitatively. The comparison was made by adapting the qualitative and quantitative methods of Sithole and Vosselman [20] for archaeology.

First, we manually classified a point cloud using context knowledge (AT, SI1, and SI2) and aerial photographs. The points were classified into the ground (2) and non-ground (1) using the Global Mapper software package, and only after careful inspection were the classifications accepted.

The qualitative assessment method is defined as a visual examination and comparison of the datasets. Five categories of circumstances in which the algorithms are likely to fail are observed:

• outliers (low or high, occur mainly due to sensor errors);
• object complexity (refers to buildings that are difficult to detect due to their size, low height, or complex shape);
• detached objects (buildings on slopes, bridges, ramps);
• vegetation (problematic categories are vegetation on slopes and low vegetation);
• discontinuity (sharp changes in a slope like cliffs and sharp ridges are treated as buildings).
  For the purposes of this study, we have added the archaeology-specific category:
• archaeological features (type 1, type 2, and type 3 as defined above in Section 2.1).

Each observed category is rated as good (item is filtered most of the time), fair (item is not filtered a few times), or poor (item is not filtered most of the time). Influence on neighboring points is also rated as good (no influence), fair (small influence), or poor (large influence).

In a similar vein, for the purposes of this study, some of the categories are irrelevant either because all filters perform good (high outliers) or because they are irrelevant to archaeology (object complexity, detached objects). Nevertheless, where present, these categories were still evaluated but are marked in grey.

It should be noted that this method is subjective by design. In this particular case, it has been carried out by two archaeologists who have extensive experience with airborne LiDAR data. It is therefore to be expected that the visual assessment is biased towards the qualities that are archaeology-specific.

The quantitative assessment was done by generating cross-matrices of each generated classified point cloud with the manually classified point cloud. The cross-matrices were then used to evaluate Type I (rejection of bare-Earth points, i.e., false negative) and Type II (acceptance of object points as bare-Earth, i.e., false positive) errors. An example of type I error is when a ground point is classified as non-ground. An example of type II error is when a non-ground point, for example belonging to a building, is classified as ground.

Quantitative assessment was carried out on small areas of 50 × 50 m within archaeological sites. If we were to apply the original method directly, the archaeology-specific accuracy would be masked by the overall accuracy (similar to how low noise was masked in the original study).

The areas were selected according to two criteria:

• the archaeological features must cover at least 50% of the area, and
• the area represents the most difficult rather than average situation.

Due to this relatively high T1 (false negative) and T2 (false positive) errors occurred in comparison to non-archaeology-specific studies, e.g., [20,47].

3. Results

3.1. Qualitative Assessment

The results of the qualitative assessment are summarized in Figure 5.

The PTIN achieved the highest average, followed by the BMHF, and PMF and SBF. The most robust performers—fewest one—were PTIN and SBF, followed by SMRF and BMHF.

On average, filters performed better on dense data than on medium and low density data. By far the best on dense data was PTIN, followed by MCC and SBF, WLS, SegBF, and BMHF. PTIN and BMHF performed best on medium density data, followed by CSF and SMRF. For low density data, PMF, PTIN, and WLS were the best.

For vegetation removal, MCC was best, followed by PMF and then WLS and BMHF as a distant third. In the preservation of discontinuities, PTIN was perfect, while SBF and WLS were acceptable. In the retention of archaeological features, SMRF and CSF were almost perfect, closely followed by SBF, PTIN, and BMHF.

PTIN was the only filter consistently among the three best; the only comparable was BMHF, which struggled with the preservation of discontinuities and failed almost completely for low-density data.

3.2. Quantitative Assessment

The results of the qualitative assessment are summarized in Figure 6. To facilitate comparison, the performance of each filter in each category was evaluated as follows. The filter with the best performance and the lowest error received nine points, the second-best eight points... the filter with the highest error—the last of nine filters evaluated—received one point. The best score in total error was obtained by SMRF, followed by a close group of BMHF, PTIN, SBF, and CSF. At the bottom of the group were PMF, SegBF, and, at some distance, WLS. The most robust performers were PTIN and BMHF (twice and once second best, respectively; never among the two worst), followed by SBF and CSF (never among the worst). The T1 error results are almost identical to total error.

The T2 error is very different, though. MCC and WLS are decidedly the best, followed by a narrow group of PTIN, BMHF, and PMF. SBF, CSF, SegBF, and at some distance SMRF are the worst.

4. Discussion

The qualitative assessment revealed that some of the categories are irrelevant for this study either because all filters perform good and would only obscure relevant differences (high outliers) or because they are irrelevant to archaeology (detached objects) or are not the subject of this article (object complexity; see archaeological features type 4 in Section 2.1 above). Nevertheless, the results are included in the table where possible (there were no detached objects in our data) to comply with the original method but are omitted from the numerical average.

Extreme examples of very dense vegetation and low vegetation (SI2) remain unsolvable for all filters. Also, the removal of vegetation on steep slopes remains problematic, but with dense data (AT) it is solved well by all filters, as is the discontinuity retention. There is a direct relationship between vegetation on the one hand and archaeology and discontinuity retention on the other hand: the more the vegetation is removed the more of archaeology and discontinuities will be removed as well.

With a few exceptions, all filters have performed well in the retention of archaeological features of types 1 (features embedded in the ground) and 2 (features partially embedded in the ground). However, there are significant differences among the filters in the costs of this retention on the vegetation removal and discontinuity retention. It is not surprising that none of the filters has fully preserved the type 3 features (standing features), as these are by definition off-terrain objects.

The results of the quantitative assessment are not clear-cut. Total error and the T1 error are almost identical, since the T1 error makes up the majority of the total error. The T2 error results are very different and—in the case of MCC, WLS, and SMRF—inverse. The latter can be easily explained: the more ground points classified (lower T1), the more are misclassified (higher T2) and vice versa. This means that MCC and WLS minimize the T2 error at the expense of the high T1 error and SMRF does the opposite. This insight could be exploited in the form of a multi-step hybrid algorithm combining, for example, MSS and SMRF. On their own MCC and WLS are suitable for high and possibly medium-density data, but not for low-density data.

On balance, BMHF and PTIN are notably better than the rest, with SBF and PMF providing acceptable performance.

5. Conclusions

Qualitative and quantitative results are not directly comparable. Thus, no filter is directly better than the rest, but the commercial filters PTIN and BMHF are consistently at the top, followed by an open access SBF. The ability of the MCC to remove vegetation is very good, as is the retention of discontinuities and archaeological features retention of the SBF. Unfortunately the BMHF—which is a hybrid of MCC and SBF—does not succeed in combining the strengths of both.

It is noticeable that in our assessment, filters based on newer algorithms (PMF, CSF, SMRF) are outperformed by some of the oldest (PTIN, SBF), which contradicts recently published findings [46]. This is a reminder of the indirect path from publishing an algorithm to the filter implementation. The latter needs to be refined and optimized, which seems to be best handled by commercial software developers. This insight puts into perspective the purpose of publishing more and more new algorithms that are only compared to the published results, e.g., [47], but are neither implemented as off-the-shelf filters nor compared to them.

Based on the qualitative assessment the archaeology-specific comparison of existing filters can be determined. The qualitative results are directly transferable to real-world problem-solving processes concerning which filter work best for a given combination of data density, landscape type, and type of archaeological features.

Concerning the influence of point density on filter performance, our results are consistent with the results presented by Sithole and Vosselman [20], namely in that the lower the point density, the worse the performance of all filters. This makes processing low-density data much more demanding.

As far as archaeological features are concerned, all filters handle the features embedded in the ground perfectly. Some filters (PMF, WLS, SegBF) struggle with semi-embedded features and all but two (SMRF, SegBF) erase standing features completely. An automatic solution for the handling of standing features has yet to be found. Therefore, a considerable manual classification effort is still inevitable [13,15].