Adaptive Framework for the Delineation of Homogeneous Forest Areas Based on LiDAR Points

Abstract

We propose a flexible framework for automated forest patch delineations that exploits a set of canopy structure features computed from airborne laser scanning (ALS) point clouds. The approach is based on an iterative subdivision of the point cloud using k-means clustering followed by an iterative merging step to tackle oversegmentation. The framework can be adapted for different applications by selecting relevant input features that best measure the intended homogeneity. In our study, the performance of the segmentation framework was tested for the delineation of forest patches with a homogeneous canopy height structure on the one hand and with similar water cycle conditions on the other. For the latter delineation, canopy components that impact interception and evapotranspiration were used, and the delineation was mainly driven by leaf area, tree functional type, and foliage density. The framework was further tested on two scenes covering a variety of forest conditions and topographies. We demonstrate that the delineated patches capture well the spatial distributions of relevant canopy features that are used for defining the homogeneity. The consistencies range from $R^2=0.84$ to $R^2=0.86$ and from $R^2=0.80$ to $R^2=0.91$ for the most relevant features in the delineation of patches with similar height structure and water cycle conditions, respectively.

1. Introduction

Forests provide a variety of ecosystem goods and services, ranging from timber production to climate mitigation [1,2]. On small scales, the canopy structure impacts the microclimate through heat fluxes and the modification of evapotranspiration and interception [1,3,4]. On a larger scale, the canopy structure also determines species richness and biodiversity, as these depend on existing niches within a forest composition [5,6,7].

In the last two decades, airborne laser scanning (ALS) has proven to be a feasible means for the derivation of forest metrics (e.g., [8,9,10]). ALS provides three-dimensional information of the vegetation distribution with high spatial resolution and with high accuracies [11,12]. Its advantage lies in the ability of the laser beam to penetrate the canopy layer [13,14,15], which enables the retrieval of detailed vertical canopy structure descriptions [16,17].

For ecological studies that investigate the above-mentioned habitat properties, a preliminary delineation of forest areas with similar conditions can be of use. Such forest areas are comparable to forest stands that form the base unit for retrieving forest inventory parameters from ALS data [8,18,19] using area-based approaches (ABA, (cf. [20,21,22])). Forest stands are thereby defined as geographically contiguous entities that are homogeneous in site type and growing stock (cf. [23]).

However, available forest stand polygons often do not reflect the actual forest characteristics but are instead based on historic forest management activities [24,25]. In our study, we therefore present a framework that enables the delineation of homogeneous forest areas from ALS point clouds. Thus, we do not use the term ‘forest stand’ in its strict meaning which, in forestry, focuses on stock- and biomass-related characteristics (cf. definition by Koivuniemi and Korhonen [23]). In contrast, we intend to widen the definition of a stand such that it is homogeneous for an application-specific characteristic. As a result of our delineation, therefore, we refer to homogeneous forest ‘patches’ rather than ‘stands’.

For delineating patches with homogeneous properties with respect to ecological process functions, such as similar canopy interception probabilities and water balance conditions, the vertical stand structure becomes explicitly relevant [4,26], and its incorporation into the delineation process is a requirement. This represents a novelty compared with existing approaches since, traditionally, forest stand delineation is performed through visual interpretation of high-resolution aerial photographs by human operators considering canopy height, canopy structure, and tree species [18,27]. Manual image analysis, however, is a time-consuming task, limited by artifacts in data as, e.g., shadow effects, prone to errors and partly subjective [18,28,29]. Previous studies, therefore, have aimed to automate forest and forest stand delineations using remote sensing data, specifically from ALS. Wang et al. [30] and Wang et al. [31] presented an approach to delineate forested from non-forested areas by the combination of ALS data and aerial images, and Straub et al. [32] and Eysn et al. [28] produced forest area delineations from ALS data that closely reflected national or international forest definitions.

Stand delineation in computer vision and image analysis amounts to a segmentation problem [27], the classic interpretation of which is the division of an image into homogeneous and spatially contiguous regions [33]. An early forest stand delineation approach using ALS data by Diedershagen et al. [18] was based on the canopy heights derived from terrain-normalized ALS point clouds. Likewise, Mustonen et al. [34] used heights from ALS-derived canopy height models (CHM) in a region-merging approach to form homogeneous stands, but they also discussed combining ALS and aerial imagery data for this purpose. Koch et al. [24] computed a raster-based classification that relies on the dominant tree species, canopy height, and variations in canopy height structure derived from ALS data. Forest stands were formed by grouping neighboring cells within each class. Sullivan et al. [35], in their approach, computed canopy cover, stem density, and average height from ALS point clouds and then used the set of these three features for an initial segmentation of homogeneous objects using a region growing approach. Because ground measurements for the stands were available, the segmented objects were then grouped into stands using supervised classification. Wu et al. [29] proposed a two-stage approach: first, the mean shift algorithm was applied to a previously computed three-band image containing bands for canopy height, density, and species, respectively, as derived from an ALS point cloud. The second stage comprised the refinement of this segmentation using regions from the region growing method computed from the three-band image.

As is apparent, these studies mainly focused on tree heights, density, maturity, and dominant tree species or types [24,27]. However, forest management concepts that favor a more diverse forest structure or succession processes on a landscape scale may lead to diversification in tree age and height, a more distinct layering, and the presence of multiple tree species within forest stands [19,36]. Delineation approaches relying predominately on height metrics will likely fail for diversified stand types where previously distinct types coexist [19]. Furthermore, our proposed method of delineating forest patches with similar conditions in ecosystem processes requires the consideration of vertical structural properties since they determine many ecosystem functions and services [37,38,39]. The feature sets used in the above-mentioned approaches, however, are close to the traditional definition of forest stands [24,29,35]: They assume a distinct and definable set of forest classes [24] or require data from field observations [35]. Although some of the frameworks potentially allow the introduction of additional information layers that describe vertical structure (e.g., Wu et al. [29]), all of these previously implemented frameworks focus on managed forests, where stands are managed in such a way that one can assume that stands are of similar age and structure, with one dominant tree type.

To meet the additional requirements and to ensure high flexibility, the homogeneity criterion in our framework is user-definable. Furthermore, no prior knowledge of the forest’s characteristics is required; the delineation approach is entirely data-driven. Moreover, it is adaptable to a range of tasks, as the attributes for forest characterization that are used in the delineation are not predefined and can be selected for specific applications.

To demonstrate the functioning of our framework, we address the following three research objectives. We

• test the framework to capture the structures present in a forest scene accurately;
• evaluate the delineated forest patches for their specificity;
• discuss the process to parameterize the algorithm.

The scope of our study is the demonstration of the ability of the framework to provide an initial set of patches with high specificity, that is, to generate a set of segments that have low variability within and large differences between segments [29,34], thus meeting the homogeneity criterion of computer vision, as specified by Cheng et al. [33]. To reveal the potential of our approach, forest patch delineation is demonstrated for two tasks. Namely, we show its applicability to the delineation of patches with similar canopy height structures. This delineation partially conforms to the strict definition of forest stands. Departing from the traditional definition, we also demonstrate the delineation of forest patches which have similar water cycle conditions. In order to confirm the robustness of the approach, the delineation was performed for test sites that cover a range of forest characteristics, including forest type, management, and species composition, and have different terrain conditions regarding topography and terrain elevation. The initial stand delineations in feature space are often spatially non-contiguous and require post-processing steps (e.g., [27,40]) in order to meet the requirement of spatial contiguity of segments in terms of image analysis [33] or stands in forestry [23,34], respectively. Our framework delivers the appropriate input for such further processing.

2. Study Data

For the study, subsets of ALS point clouds from country-wide acquisitions were selected, covering a variety of forest and site conditions.

2.1. Burgenland Scene

The Burgenland scene covers a 1200 × 2400 m area (northing × easting) centered at 5280900 N, 606300 E (UTM33 N; Figure 1). The scene covers terrain heights between 415 and 590 m a.s.l. The 95%-quantile of vegetation height is 24.1 m. As the canopy height map in Figure 1 illustrates, the scene comprises a wide range of canopy heights, including clear-cut areas. The variation in height can also be discerned from the canopy profiles in Figure 1 that further illustrate the large variation in canopy thickness. Furthermore, the scene covers a wide diversity regarding canopy cover, canopy density, and tree functional types. ALS data were acquired in April 2010 under leaf-off and snow-free conditions using Riegl LMS-Q560 and LMS-Q680 sensors. The mean flying height was 400 m above ground. The data were acquired with a median pulse density of 15 pulses/m${}^2$ from which the maximum number of returns per pulse is 9. The resulting median point density derived from the data is 24 points/m${}^2$.

2.2. Ötscher Scene

The Ötscher scene covers a 1415 × 1550 m area centered at 304540 N,-84970 E (MGI/Austria GK East; Figure 2) and is characterized by a large elevation gradient with terrain heights ranging from 800 to 1650 m a.s.l., including the upper timberline in the southern parts of the scene. The 95%-quantile of vegetation height is 28.8 m. The study area covers a variety of forest types regarding tree heights, canopy density, and canopy thickness. ALS data were acquired in January 2007 in leaf-off conditions using a Riegl LMS-Q560 sensor with a mean flying height of 620 m above ground. The data were acquired with a median pulse density of 12 pulses/m${}^2$. The maximum number of returns per pulse is 6. The resulting median point density within the scene is 19 points/m${}^2$.

3. Methods

As depicted in Figure 3, the framework developed in this study comprises an initial feature computation step (denoted by Feature calculation in Figure 3), followed by the actual patch delineation step (denoted by Iterative Splitting Segmentation framework in Figure 3). The patch delineation framework is a closed pipeline and consists of three major steps, namely, (i) a splitting step, (ii) an elimination step in which clusters that are too small are eliminated, and (iii) a final merging step in which overlapping clusters are merged to a final forest patch class. Each of these three steps is iteratively run, and the next processing step is initiated if a criterion is fulfilled. The entire delineation is performed in feature space in which, for the control of each of the three steps (indicated by Evaluation in the figure), a different feature (indicated by Feature in the figure) and different thresholds (indicated by Threshold) can be defined. The result is a labeled point cloud, where a label specifies the forest patch class to which the respective point belongs.

3.1. Feature Computation

3.1.1. Forest Structure Characterization

We computed seven point-cloud features from the terrain-normalized point cloud. Features were computed from all points (i.e., LiDAR echoes) lying within a cylindrical neighborhood centered at a query point from which radii of 2, 5, and 10 m were used for the neighborhood search. In the following, the features describing local canopy characteristics are discussed in detail.

Fractional Cover

As discussed in Morsdorf et al. [14], fractional cover is a dimensionless parameter describing the fraction of ground that is covered by vegetation within a certain ground area, and it can be used to measure crown coverage within a point’s neighborhood. The ratio is computed from the point cloud as:

$$fractionalcover=\frac{N_{canopy}}{N_{total}}$$

(1)

where $N_{canopy}$ denotes the number of points derived from the canopy, and $N_{total}$ is the total number of points in the respective neighborhood, i.e., canopy and non-canopy echoes. In this study, a canopy threshold of 3 m above terrain was used to distinguish canopy from non-canopy points.

Canopy Density

The canopy densities d40 and d50 were computed following the method in Hollaus et al. [41], as follows:

$$dQQ=\frac{N_{hQQ}}{N_{total}}$$

(2)

where $N_{hQQ}$ corresponds to the number of points above a certain height quantile with regard to the maximum point height within a neighborhood cylinder ($40\%$- and $50\%$-quantile, respectively, in our case), and $N_{total}$ corresponds to the total number of points within that neighborhood. Hollaus et al. [41] found the d50-feature to be well suited to distinguishing coniferous from deciduous trees.

95%-Height Quantile

In order to have a height attribute at hand that is robust to outliers, we computed canopy height as the 95%-point height (h95) within a given neighborhood.

Vegetation Profiles

Vegetation profiles from discrete ALS data were derived following the approach by Coops et al. [17]. The computation consists of first calculating a gap probability profile $P_{gap}\left(h\right)$:

$$P_{gap}\left(h\right)=1-\frac{N_{\{h>z\}}}{N_{total}}.$$

(3)

$N_{\{h>z\}}$ denotes the number of points above height z, where a bin height $dz$ of 1 m was used as the step size to iterate over the entire canopy height. The gap probability profile is then transformed into the cumulative projected plant area index by $PAI\left(h\right)=-ln\left(P_{gap}\left(h\right)\right)$, indicating the plant area index (PAI) above height z above ground. The subsequent formation of the first derivative of the PAI profile results in the vegetation profile. Following the method by Palace et al. [42], from this profile, we interpret the number of layers (vg_nLayers) as the number of peaks, the height of the highest peak as the topmost layer (vg_topLayer), and the distance between the top and the bottom layer peaks (vg_layerDiff) as the canopy thickness.

3.1.2. Introduction of Scale

Structure descriptors are required to characterize the forest structure on an appropriate scale. For the computation of features describing the local forest structure around a point, all points lying within a certain point’s neighborhood are used. Contrary to urban applications for which the optimum neighborhood can be derived from the point cloud [43], the scale in forestry applications has to match the extent of the objects, i.e., the trees. For the features utilized in this study, search radii of 2–15 m [14] or squared neighborhoods with side lengths between 1 and 10 m [9,17,44] were deployed as in previous studies.

Segmentation of forest stands requires the generalization of the fine-scale variabilities that are introduced by, for example, small gaps between trees. Diedershagen et al. [18] already discussed how such small-scale heterogeneities can hamper proper stand delineations. Therefore, the aim is to regard such small-scale heterogeneities as uniformly homogeneous if they are similar across a larger area. We, thus, chose a search radius of 10 m around each point as the maximum neighborhood size, which exceeds the scale of individual tree crowns and lowers the impact of small-scale gaps. Along with this large neighborhood, we also computed the features for search radii of 2 and 5 m in order to retain a certain amount of variance at smaller scales. The computation of the described neighborhood search within the point cloud was accomplished using the kd-tree search implemented in scipy.spatial from the SciPy library [45].

The drawback of large search radii is an exhaustive neighborhood search within the point cloud. To reduce the computation cost, we did not calculate features for every LiDAR echo contained in the point cloud but performed a preliminary sampling of query points in a regular grid with a 1 m side length using point cloud tools by Glira [46]. Structure features were only computed for the LiDAR echo (an actual element of the ALS point cloud) closest to the cell center for each grid cell while exploiting the entire, non-sampled point cloud. The delineation is subsequently performed on the sampled point cloud, in which each point has a number of features. The envisaged patch delineation is a map with cell sizes that are a multiple of this sampling interval. We used this procedure to avoid reducing the resolution in a way that would affect the final result.

3.2. Iterative Splitting Segmentation

Our framework follows the divide-and-conquer paradigm [47], yet it is only in partial agreement with the original concept of this class of algorithms. The point cloud is iteratively split into two subsegments. Then, independently for each subsegment, a decision is made whether to continue the splitting or stop the processing. In contrast to the original algorithm paradigm, however, these subsegments are not subproblems that need to be recombined by necessity. That is, subsegments are only recombined if the size of a subsegment undercuts a certain threshold or if a subsegment is highly similar to another subsegment.

The advantage of this approach is that no preliminary knowledge of the number of final forest patch classes is required, instead the number is found from the data based on a user-defined homogeneity criterion.

3.2.1. Splitting Step

The splitting step performs an iterative bipartitioning of a set of points into two subsets using the k-means algorithm [48] implemented in MATLAB 2017b [49]. Initially, the entire point cloud is employed for k-means clustering. Since k-means clustering uses Euclidean distance to assign points to cluster centers, all attributes were scaled to have the range $[0,1]$.

The idea is to iteratively bipartition the point cloud as long as two criteria are fulfilled. First, the splitting of the point set Sub0 is only continued if it contains a minimum number of points defined by the threshold $th\_minN\_split$. If this criterion is met, Sub0 is passed to the k-means clustering step and the resulting two subsets $\left\{Sub1,Sub2\right\}$ are evaluated for their uniqueness based on a user-defined attribute $ftr\_thSplit$. The overlap of attribute $ftr\_thSplit$ in the feature space of the two subsets is evaluated using the Mahalanobis distance metrics:

$keepSplit(Sub1,Sub2)=\{\begin{array}{ll}1,& if[mdis{t}_{Sub2}\left(mean\left(Sub1\right)\right)>th\_split]\&\\ & [mdis{t}_{Sub1}\left(mean\left(Sub2\right)\right)>th\_split]\\ 0,& \mathrm{otherwise}\end{array}$ where $mdis{t}_{Sub2}\left(mean\left(Sub1\right)\right)$ measures the Mahalanobis distance of the mean of Sub1 from the distribution of Sub2, which is the distance in multiples of standard deviations of the distribution of Sub2 [50]. If the Mahalanobis distances exceed the threshold $th\_split$, we consider $\left\{Sub1,Sub2\right\}$ to be unique in the relevant feature $ftr\_thSplit$, and the two subsets re-enter the splitting stage independently. If the condition is violated, the performed subdivision is undone and further subdivision of Sub0 is terminated.

Typically, the splitting step results in an oversegmentation of the initial point cloud into a too large a number of clusters. The implication is that the resulting clusters are potentially too small regarding their number of contained points and regarding the cluster uniqueness in feature space. Therefore, two further steps are necessary to ensure an adequate generalization of the patch delineation and to reduce the number of homogeneous classes to a level that is appropriate for a respective application.

3.2.2. Elimination of Small Clusters

The first generalization step consists of the treatment of clusters that are too small (i.e., the point sets from the splitting step), where the threshold $th\_minSegSize$ defines the minimum number of points to be contained in a final stand class. This implies that too-small clusters need to be merged with other clusters, which is again done in the feature space. The approach is to join a cluster that is too small to the cluster for which the change in homogeneity after merging is the smallest in a user-defined attribute $ftr\_mergeSegments$ (Equation (4)).

$$\underset{i}{argmin}|std\left(Se{t}_i\right)-\left(std(Se{t}_{small}\cap Se{t}_i)\right)|$$

(4)

Starting with the smallest cluster to be merged, clusters being too small $Se{t}_{small}$ are iteratively compared with all other clusters $Se{t}_i$ from the splitting step, including other clusters containing too few points. The merging step is terminated as soon as the number of points contained in each cluster exceeds $th\_minSegSize$.

3.2.3. Elimination of Non-Unique Clusters

Although clusters comprise a minimum number of points, their uniqueness is not guaranteed per se, and clusters may considerably overlap in the feature space. In a final step, the cluster overlap in a user-defined attribute $ftr\_mergeOverlap$ is evaluated, which is again measured as the Mahalanobis distance of $mean\left(Set1\right)$ to the distribution of $Set2$ and vice versa. If the overlap of two clusters in both directions falls below a threshold $th\_mergeOverlap$, the two respective clusters are merged.

3.3. Validation

The validation of the segmentation framework is threefold in order to address the stated research objectives. First, the delineation framework was tested for two specific applications in order to demonstrate that the framework appropriately captures and reflects the canopy structural properties. Second, different levels of segmentation were compared regarding the obtained specificity of the patches, and third, the delineation consistency and its dependence on the applied thresholds was investigated using sensitivity analysis in order to propose a process for their parameterization.

3.3.1. Forest Patch Delineation

The segmentation task requires a set of point cloud features that capture forest characteristics. The final forest classes correspond to homogeneous clusters within that feature space. As no ground truth data for forest characteristics was available for the two scenes, the segmentation framework was first assessed by visually comparing the spatial distribution of the patches with the spatial distributions of attributes relevant to the delineation under study. This indicates whether the segmentation framework is able to capture the relevant forest characteristics and represent breaks between patch classes. Second, attribute distributions within and between forest patch classes were assessed (i) by comparing the value distributions of relevant attributes and (ii) by statistical analysis of the segmentation consistencies using the $R^2$ metric. This metric has been used by Mustonen et al. [34] and Wu et al. [29] in stand delineation evaluations to compare the delineation results with ground-measured canopy metrics. In our study, we used it to evaluate within-class homogeneities, so the measure is turned into a consistency rather than an accuracy measure. $R^2$ was computed as

$$R^2=1-\frac{S{S}_{within}}{S{S}_{total}}.$$

(5)

where $S{S}_{within}$ captures the variability within forest classes as

$$S{S}_{within}=\sum _{i=1}^k\sum _{j=1}^{n_i}{(x_{ij}-\overline{x_i})}^2,$$

(6)

and $S{S}_{total}$ measures the total variability within the scene as

$$S{S}_{total}=\sum _{i=1}^k\sum _{j=1}^{n_i}{(x_{ij}-\overline{x})}^2$$

(7)

with k denoting the number of forest classes, $n_i$ the total number of points within forest class i, $x_{ij}$ the attribute value of point j in forest class i, $\overline{x_i}$ the mean attribute value of all points within class i, and $\overline{x}$ the mean attribute value of all points within the scene. $R^2$ has a range of [0 1], and interpretation is such that the higher the $R^2$ value, the more consistent is the segmentation [29].

The performance of the segmentation framework was tested for the segmentation of forest patches with a similar height structure (Experiment 1) and for the delineation of forest patches with homogeneous water cycle conditions (Experiment 2). Experiment 1 additionally allowed for visual inspection of forest profiles, so we analyzed whether differences in canopy structure are recognizable between the different patch classes.

Experiment 1: Forest Patch Segmentation with Similar Height Structure

In Experiment 1, we selected a set of features related to the height structure and canopy thickness of the forest. The fulfillment of the homogeneity criterion in these canopy characteristics is well recognizable also from visual inspection of forest profiles and allows for visual evaluation in addition to statistical metrics alone.

Tree height alone does not cover all facets that are commonly used for stand delineations. Previous studies aiming to delineate stands with similar growing stock also exploited species and crown coverage for the stand characterization [24,29,35], and such properties are missed if analyzing the canopy height alone. Yet, studies by Diedershagen et al. [18] and Mustonen et al. [34] revealed the consistency between the results of delineation approaches relying on height metrics alone and those of manual stand delineations.

From the computed feature set, we selected the h95, $vg\_layerDiff$, $vg\_topLayer$, and $vg\_nLayers$ attributes for k-means clustering. Furthermore, different scales of neighborhoods were considered for their computation (

Table 1. Point cloud features used in the segmentation framework for Experiment 1.

Usage in Segmentation Pipeline	Point Cloud Feature
Input to k-means	h95 for search radii = {2, 5, 10 m}
	$vg\_layerDiff$ for search radii = {2, 5, 10 m}
	$vg\_topLayer$ for search radii = {2, 5, 10 m}
	$vg\_nLayers$ for search radius = {10 m}
$ftr\_thSplit$	h95, search radius = 5 m
$ftr\_mergeSegments$	h95, search radius = 5 m
$ftr\_mergeOverlap$	h95, search radius = 5 m

). As we considered canopy height to be the most important property in our approach to delineating patches, we chose the h95-attribute as a criterion for the threshold evaluations.

Experiment 2: Forest Patch Segmentation for Water Cycle Studies

The canopy affects the water cycle through interception and evapotranspiration. The interception capacity of a canopy is thereby a function of its leaf area, tree functional type, and foliage density [4,51] while, for canopy transpiration, leaf area is the key canopy element [4,26]. However, most of these canopy parameters cannot be derived directly from the point cloud but require a set of ALS metrics that can be used as proxies. Therefore, we translated the required canopy metrics into point cloud features, as reported in

Table 2. Canopy parameters with relevance to water cycle properties and the respective airborne laser scanning (ALS) point cloud features used as proxies.

Canopy Parameters	ALS Metrics
Leaf area	fractional cover, $vg\_nLayers$
Foliage density	d40, $vg\_layerDiff$
Tree functional type	d50

, that again were computed at different scales (

Table 3. Point cloud features used in the segmentation framework for Experiment 2.

Usage in Segmentation Pipeline	Point Cloud Feature
Input to k-means	fractional cover for search radii = {2, 5, 10 m}
	d40 for search radii = {2, 5, 10 m}
	d50 for search radii = {2, 5, 10 m}
	$vg\_layerDiff$ for search radius = {2, 5, 10 m}
	$vg\_nLayers$ for search radius = {10 m}
$ftr\_thSplit$	fractional cover, search radius = 5 m
$ftr\_mergeSegments$	d50, search radius = 5 m
$ftr\_mergeOverlap$	d50, search radius = 5 m

). In this experiment, metrics related to leaf area (fractional cover) and to species (d50) were employed to control the segmentation pipeline.

3.3.2. Definition of Number of Patch Classes

A question that arises in segmentation frameworks is regarding the appropriate number of resulting patch classes into which the point cloud ought to be delineated. To approach that aspect, the adequacy of results was evaluated (i) through the analysis of the specificity of patch class distributions of relevant features, (ii) through the segmentation consistency metrics $R^2$, and (iii) by inspection of the sizes of the resulting forest patches.

3.3.3. Sensitivity to Thresholds

The definitions of $th\_mergeOverlap$ and $th\_split$ have a large impact on the final delineation results. To test the sensitivity of the delineation to the selection of the respective threshold values, the delineation results for a number of parameterizations were tested and compared.

4. Results

4.1. Feature Computation

The computation of features for a subsample of the initial points with a sampling interval of 1 m reduced the number of points to $3.8\%$ and $11.3\%$ of the initial point cloud size for the Burgenland and Ötscher scenes, respectively. Maps of the spatial distribution for three features fractional cover, return height ratio d40, and canopy thickness vg_layerDiff are depicted in Figure 4 and Figure 5 for the Burgenland and Ötscher scene, respectively. Figure 4 further illustrates the impact of the search radius on the attribute computation.

4.2. Iterative Splitting Segmentation

The functioning of the iterative splitting step is illustrated in Figure 6. Starting with the unsegmented point cloud, the maps in the top row depict the spatial distribution of the two resulting segments after each bipartitioning step for iterations 1–3. For a better illustration, only the splitting of the larger of the two clusters (colored in blue in each map) in terms of the respective point count is traced, while the smaller cluster (colored in red) after each splitting step is excluded from the visualization in the subsequent iteration. The scatterplots in Figure 6 (bottom row) show the attribute distributions of the two clusters in the respective maps in the h95–vg_layerDiff feature space after each splitting iteration.

4.3. Validation

4.3.1. Forest Patch Delineation

Experiment 1

The forest patch delineation for areas with similar height structure (Experiment 1) uses the features listed in Table 1. Figure 7 and Figure 8 depict the resulting forest delineations employing the parameterizations specified in

Table 4. Parameterization of thresholds used for the feature set of Experiment 1 for the Burgenland scene with $th\_split$ held constant while varying $th\_mergeOverlap$. Further stated is the consistency $R^2$ of the resulting delineations for the canopy height (h95) and canopy thickness ($vg\_layerDiff$), whose homogeneity we considered the most important for the task.

Threshold	Number of Segments k
	k =4	k =5	k =8
$th\_split$	0.4
$th\_minN\_Split$	1000
$th\_minSegSize$	20000
$th\_mergeOverlap$	5.0	3.0	1.3
$R_{h95}^2$	0.79	0.84	0.85
$R_{vg\_layerDiff}^2$	0.70	0.77	0.78

and

Table 5. Parameterization of thresholds used for the feature set of Experiment 1 for the Ötscher scene with $th\_split$ held constant while varying $th\_mergeOverlap$. Further stated is the consistency $R^2$ of the resulting delineations for the canopy height (h95) and canopy thickness ($vg\_layerDiff$), whose homogeneity we considered the most important for the task.

Threshold	Number of Segments k
	k =4	k = 5	k = 7
$th\_split$	0.4
$th\_minN\_Split$	1000
$th\_minSegSize$	20000
$th\_mergeOverlap$	2.3	2.0	1.5
$R_{h95}^2$	0.75	0.86	0.86
$R_{vg\_layerDiff}^2$	0.60	0.68	0.68

for the Burgenland and Ötscher scenes, respectively. Depicted are the maps for the delineations resulting in $k=5$ and $k=4$ patch classes for both scenes. The attribute value distributions of the resulting patches are shown in Figure 9 (Burgenland) and Figure 10 (Ötscher). Additionally, forest profiles from the two scenes for the segmentation into $k=5$ patches are displayed in Figure 11 (Burgenland) and Figure 12 (Ötscher). The consistencies $R^2$ of the patch delineations are listed in Table 4 (Burgenland) and Table 5 (Ötscher).

Experiment 2

For the delineation of forest patches with similar water cycle conditions, the segmentation framework was fed the feature set reported in Table 3 and parameterized as stated in

Table 6. Parameterization of thresholds used for the feature set of Experiment 2 in the Burgenland scene. Further stated is the consistency $R^2$ of the resulting delineations for fractional cover ($fc$) and canopy density (d50), whose homogeneity we considered the most important for the task.

Threshold	Number of Segments k
	k =4	k = 5	k =6
$th\_split$	0.22
$th\_minN\_Split$	1000
$th\_minSegSize$	20,000
$th\_mergeOverlap$	1.2	0.95	0.8
$R_{fc}^2$	0.67	0.80	0.82
$R_{d50}^2$	0.65	0.76	0.79

and

Table 7. Parameterization of thresholds used for the feature set of Experiment 2 in the Ötscher scene. Further reported is the consistency $R^2$ of the resulting delineations for fractional cover ($fc$) and canopy density (d50), whose homogeneity we considered the most important for the task.

Threshold	Number of Segments k
	k = 4	k = 5	k = 6
$th\_split$	0.4
$th\_minN\_Split$	1000
$th\_minSegSize$	20,000
$th\_mergeOverlap$	4.0	2.5	0.7
$R_{fc}^2$	0.90	0.91	0.91
$R_{d50}^2$	0.88	0.90	0.90

for the Burgenland and Ötscher scenes, respectively. Figure 13 and Figure 14 show the respective spatial distributions of the homogeneous patches for the two scenes, while the attribute distributions of the classes are depicted in the boxplots in Figure 15 and Figure 16, respectively.

4.3.2. Definition of Number of Patch Classes

Table 8. Resulting segment sizes for segmentation into $k=4$, $k=5$, and $k=8$ patch classes for the Burgenland scene, Experiment 1. The sampled point cloud contained 2,874,780 points for which canopy features were calculated. All values are in [%] with regard to this point cloud size. The color scheme follows that in the maps in Figure 7 while the three additional classes S3, S7, S8 are derived from segmentation into $k=8$ classes.

Class Label	Number of Classes k
	k = 4	k = 5	k = 8
blue/S1	3.44	3.44	3.44
orange/S2	41.12	10.55	7.14
purple/S5	34.71	34.71	34.71
green/S6	20.73	20.73	18.24
yellow/S4		30.57	30.57
S3			3.42
S7			0.70
S8			1.79

and

Table 9. Resulting segment sizes for segmentation into $k=4$, $k=5$, and $k=7$ stand classes for the Ötscher scene, Experiment 1. The sampled point cloud contained 1,983,117 points for which the canopy features were calculated. All values are in [%] with regard to this point cloud size. The color scheme follows the maps in Figure 8 while the two additional classes, labeled orange and light blue, are derived from segmentation into $k=7$ classes.

Class Label	Number of Classes k
	k = 4	k = 5	k = 7
blue	8.15	8.15	8.15
yellow	12.93	12.93	11.00
green	66.11	57.41	57.41
red	12.80	12.80	10.49
purple		8.70	8.70
orange			1.93
light blue			2.31

report the segment sizes of the final forest patch classes in relation to the total number of points for $k=4$, $k=5$, and $k=8$ classes for Burgenland and for $k=4$, $k=5$, and $k=7$ for the Ötscher scene, respectively, for Experiment 1. For the Burgenland scene, the additional solution for the $k=8$ patch classes resulted from the threshold $th\_mergeOverlap=1.3$ while, for Ötscher scene, the $k=7$ classes were achieved with a threshold of $th\_mergeOverlap=1.5$. In combination with the respective consistency measures $R^2$ in Figure 4 (Burgenland) and Figure 5 (Ötscher), this allows for the detection of an appropriate number of forest patch classes.

4.3.3. Sensitivity to Thresholds

Related to the question of the number of desired forest patch classes is the definition of the thresholds $th\_meanOverlap$ and $th\_split$. The sensitivity of the resulting number of forest patch classes to these thresholds was therefore analyzed.

Sensitivity to $th\_mergeOverlap$

The impact of $th\_mergeOverlap$ on the resulting number of patch classes was analyzed using a constant value for $th\_split$ and varying $th\_meanOverlap$ in the range 1.10–5.42.

Table 10. Resulting number of patch classes k depending on the parameterization of the merging threshold th_mergeOverlap. Depicted are the results for the Burgenland scene using the feature set from Table 1 and the thresholds depicted in Table 4.

Number of Classes	3	4	5	6	7	8	9
$th\_meanOverlap<$	5.42	5.00	4.09	2.98	1.79	1.75	1.10

reports an exemplary number of resulting patch classes depending on $th\_mergeOverlap$ for Experiment 1 for the Burgenland scene, where the number of classes prior to the merging step was $k=17$. The attributes exploited for k-means clustering and features used as thresholds are reported in Table 1. As is visible, the increase (i.e., delta) in $th\_mergeOverlap$ that is required to merge two additional classes into a single one varies in a range from $k=9$ to $k=3$ classes. While the results with $k=9$ to $k=7$ classes and $k=4$ to $k=3$ classes are close together in terms of differences in $th\_mergeOverlap$ (the difference in $th\_mergeOverlap$ is 0.69 and 0.42 between $k=9$ and $k=7$ and $k=4$ and $k=3$ classes, respectively), the required differences in $th\_mergeOverlap$ are larger in the span between the $k=7$ and $k=4$ classes.

Sensitivity to $th\_split$

To investigate the impact of $th\_split$ on the resulting segmentation, $th\_split$ was varied while all other parameters were held constant. The sensitivity analysis again was performed on the Burgenland scene using the parameterization from Experiment 1 (Table 1).

Table 11. Sensitivity of the segmentation framework to the parameter $th\_split$ for the Burgenland scene. The table reports the chosen threshold values and the resulting number of classes after the three processing steps comprising the iterative splitting step (k-means), the elimination of small clusters (elimination), and the merging of overlapping clusters (overlap merging).

Threshold	Parameter Value
$th\_split$	0.1	0.4	1.0
$th\_minN\_Split$	1000
$th\_minSegSize$	20,000
$th\_mergeOverlap$	3.0
Processing Step	Number of Classes
k-means	217	12	7
elimination	33	10	7
overlap merging	5	5	7

depicts the parameters used for the thresholds. Furthermore, the table states the resulting number of patch classes after the initial splitting step (denoted by k-means), after the elimination of too-small clusters (denoted by elimination), and after merging overlapping clusters in the feature space (denoted by overlap merging) for $th\_split$ values of 0.1, 0.4, and 1.0, respectively.

Table 12. Delineation consistency $R^2$ for the canopy height (h95) and the canopy thickness ($vg\_layerDiff$) measured as $R^2$ for different $th\_split$ values for the Burgenland scene after merging into $k=5$ classes.

Threshold	Parameter Value
$th\_split$	0.1	0.4	1.0
$th\_mereOverlap$	3.0	3.0	5.3
Attribute	Consistency${\mathit{R}}^{\mathbf{2}}$
h95	0.60	0.84	0.80
$vg\_layerDiff$	0.47	0.77	0.72

reports the consistencies in the h95 and the vg_layerDiff attributes that result if the respective delineations are processed until a class number of $k=5$ is reached. The maps in Figure 17 depict the resulting forest patch delineations for $th\_split=0.1$ (left) and for $th\_split=1.0$ (right), respectively (note that in the first map, the scene was delineated into $k=4$ classes instead of $k=5$, as the size of the fifth class was very small).

5. Discussion

5.1. Feature Computation

Feature computation is the most computationally intensive step, also because of the three different neighborhoods (2 m, 5 m, 10 m). The main effect of larger search radii is the spatial smoothing of the attribute distributions (Figure 4), leading to smoothed boundaries between the forest areas. We found the combined exploitation of these three neighborhoods to be necessary, as they contribute to segmentation in different ways. Small neighborhoods reflect small-scale heterogeneities and are required to trace sharp borders between different forest patches. Larger neighborhoods, on the other hand, are a prerequisite to ensure generalization, as the forest structure description should go beyond single tree crowns. Larger search radii, however, lead to transition zones at sharp forest type boundaries, due to averaging over larger areas in the feature calculation.

The chosen approach, in which features are computed for a regularly subsampled point cloud, provides a level of detail that includes sufficient structural details required for the task of large-scale forest delineations while still reducing the computation time. Since the actual calculation of features exploits the initial point cloud, the spatial consistency of the attribute fields can be ensured (Figure 4 and Figure 5).

Moreover, as features are calculated from the entire point cloud, the canopy penetration properties are well captured. This allows for the computation of features that have already been used in studies which focused on describing the vertical canopy structure (e.g., [14,52,53]), as well as for stand delineations [29]. Vertical structure descriptors are a prerequisite for the task of, for instance, forest patch delineations in the context of water cycle studies.

5.2. Iterative Splitting Segmentation

The aim of this study is the development of a segmentation framework on point entities, where each point is attributed a number of features describing the local canopy characteristics. The delineation of homogeneous forest patches is a complex task for several reasons: First, potentially high variability present in a forest scene should be grouped into a small number of distinct patches; second, forest habitats may be characterized by a continuous transition rather than by hard breaks, especially in case of natural forests; and third, gaps and holes constitute small-scale heterogeneities within otherwise homogeneous patches [18] and therefore should be ignored.

To ensure a coherent patch delineation without an excessive number of resulting clusters in the feature space, we propose a top-down segmentation approach. This is in contrast to bottom-up approaches used in, for example, Wang et al. [54], where the sought classes differ more distinctly for a certain property. In such region growing approaches, segmentation is performed by subsequently joining a set of homogeneous entities obtained through prior oversegmentation if the homogeneity criterion is met.

The aim of the first step in our delineation framework is the realization of a set of clusters with a large within-cluster homogeneity and a large dissimilarity between the clusters. As the scatter plots in Figure 6 (bottom row) depict, the forest areas in each cluster are increasingly more homogeneous after each splitting iteration. Conversely, this means that within-cluster differences become increasingly less pronounced after each bipartitioning of a subsegment. Therefore, the splitting iteration is finally terminated. The result of this initial step is typically a slight oversegmentation, where the degree of oversegmentation is user-specified (threshold $th\_split$).

Despite the lack of an explicit spatial constraint being introduced into the segmentation framework, the splitting step results in a set of segments which are well connected in space (Figure 6). The high spatial consistency is a result of both the segmentation concept as well as the semantic meaning of the employed feature set. The implemented top-down segmentation first splits the forest scene along the largest differences in the feature space, and minor differences are only gradually considered. As the features used for segmentation capture structural properties of the forest that are relatively continuous across space (Figure 4 and Figure 5), spatially adjacent points are more closely grouped in the feature space than more remote ones. We thus deduced that iterative bipartitioning, which represents a top-down segmentation approach, with the exploitation of structural features from different scales is a suitable framework for the initial delineation of homogeneous forest patches.

The initial splitting step, however, is likely to produce segments that are too small. Therefore, following the approach by Wang et al. [54], we performed a subsequent merging step in the feature space wherein clusters with similar structural canopy properties are merged.

Although the specified minimum segment size criterion is met after this step, the number of clusters still tends to be too large (Table 11). Furthermore, we found some of the clusters to have considerable overlap in the feature space of the relevant attributes. Thus, in order to increase the benefit of patch delineations for forestry or ecological applications, it might be reasonable to merge such overlapping clusters. We decided to deploy the Mahalanobis distance to detect cluster overlaps, as it indicates the distance of the cluster means with respect to the standard deviations of the attribute distributions. This is crucial because the goal is to not merge clusters that are close in the feature space but to have little overlap. However, in order to offer the user the greatest possible flexibility in this final step, the rigor of merging such overlapping clusters is user-specifiable by adjusting the threshold $th\_meanOverlap$, and it also allows for leaving the delineation solution unchanged after the initial merging, where only small segments are eliminated.

Finally, the point entity-based patch delineation provides the possibility to define the desired scale in a post-delineation step, where the delineated point cloud typically is transformed into a raster structure. This also contributes to the flexibility of the approach and is in contrast to previous studies that worked on raster images, where features are already computed from an initially rasterized CHM (cf. [24,35]), or that exploited the point cloud but transformed the derived attributes into a raster matrix for the actual segmentation (cf. [29]). In comparison, limitations to the spatial resolution of the delineation in our framework are only introduced by the sampling rate of the points for which features are computed. Delineation using a narrower point sampling scheme or even using the original point cloud would also be conceivable. These limits are only applied due to the computation costs of the feature calculation.

The implementation of the delineation pipeline, therefore, allows the solution to be adapted to a user-defined application by specifying the tolerated cluster overlap, the number of patch classes, and the output resolution. The high consistencies represented by the $R^2$ metrics, in our view, justify the approach.

5.3. Validation

5.3.1. Forest Patch Delineation

Contrary to previous studies that exploited a predefined set of features for patch delineations and performed tests using managed forests characterized by even-aged trees and one dominant tree species (e.g., [18,24,29]), our delineation framework is not explicitly based on a certain set of predefined features but rather is very flexible regarding the attributes that are exploited for the patch delineation. The experiments carried out in the study prove that the delineation framework is functional for two different applications. We could further demonstrate our approach for a variety of forest types and for forests under different topographic conditions.

Experiment 1

The delineation of forest areas with similar canopy height structure generates patches that well reflect the spatial variations of the h95 and $vg\_layerDiff$ attributes. The patches further reveal distinct differences in the profiles (Figure 11 and Figure 12) and in the attribute value distributions (Figure 9 and Figure 10), where the median values and the 25th –75th percentile ranges of the attributes show a distinct grading across the classes.

Analysing the data visually, the five classes of the Burgenland scene can be interpreted as dense, juvenile forest (blue), tall forest with distinct multi-layering (green), and a slightly less tall forest class with a smaller number of distinct layers (purple). The other two classes (orange and yellow) represent medium forest of different height. In the Ötscher scene, the classes are: Tallest, multi-layered forest (red), medium height (green), small-grown trees (yellow), transition between yellow and green class, with a shorter canopy (purple), and low-growing vegetation (blue), mainly in the highest regions (southern part).

The $R^2$ metrics in Table 4 and in Table 5 indicate that the consistency is higher for the delineation into $k=5$ classes ($R^2=0.84$ and $R^2=0.86$ for the Burgenland and the Ötscher scenes, respectively) than for the delineation into $k=4$ classes ($R^2=0.79$ and $R^2=0.75$ for the Burgenland and the Ötscher scenes, respectively) when the h95 attribute is considered. The consistency for the $vg\_layerDiff$ attribute in general is lower in both scenes ($R^2=0.77$ and $R^2=0.70$ for $k=5$ and $k=4$ classes, respectively, for the Burgenland scene; $R^2=0.68$ and $R^2=0.60$ for $k=5$ and $k=4$ classes, respectively, for the Ötscher scene).

Experiment 2

The spatial distribution of the patch classes with similar water cycle conditions again reflects the spatial distribution of the attributes (Figure 4 and Figure 5). The boxplots in Figure 15 and Figure 16 show good separability of the delineated patches when the median values and the value ranges of the 25th–75th percentiles are considered. The high degree of homogeneity that is achieved is reflected by the values of the consistency metric $R^2$. For the two most important attributes and the delineation into $k=5$ classes, the values are $R^2=0.80$ and $R^2=0.91$ for fractional cover for the Burgenland and Ötscher scenes, respectively, and $R^2=0.76$ and $R^2=0.90$ for d50 for the Burgenland and Ötscher scenes, respectively.

When the merging of overlapping clusters is performed more rigorously, resulting in $k=4$ patch classes, the consistency decreases to $R^2=0.67$ and $R^2=0.65$ for fractional cover and d50, respectively, for the Burgenland scene, and to $R^2=0.90$ and $R^2=0.88$ for fractional cover and d50, respectively, for the Ötscher scene. For this more rigorous merging, the two classes with the highest values for the d50 and fractional cover attributes are unified in both scenes. On the other hand, contrary to the two discussed attributes, the boxplots reveal that the canopy thickness (vg_layerDiff) is not a unique characteristic between the patch classes, as the class overlap is large (except for low vegetation, colored in blue in both scenes).

A recognized deficiency that is inherent in the chosen approach and applies to both experiments is its inability to capture sharp forest borders. Firstly, this is due to the large neighborhoods that were used for feature computation, which smooths out the attribute fields in space, making sharp structural boundaries less distinct (Figure 4, top row). This poses a problem because, secondly, the entire patch delineation is performed in the feature space without explicit consideration of either the point locations or the attribute values of spatially adjacent points. Large neighborhood definitions, however, cannot be neglected in the segmentation, as this would increase the importance of the small-scale variabilities present in the forest structure, such as those resulting from the gaps between tree crowns. This would hamper patch delineations [18]. Comparable issues to inappropriate borders were also reported for the forest delineation approach by Eysn et al. [28]. As in their approach, the correction of the border lines needs be done in the object space. In our framework, the width of the transition lines is in relation to the search radii that were used in the feature computation.

An open aspect to completely fulfill the strict patch definition by Koivuniemi and Korhonen [23] is in relation to the spatial consistency of the delineated patches. This is not yet given, and the results still contain speckle that, depending on the intended further use of the delineated patches, should be eliminated. An overview of the process to tackle this issue is given in Schindler [40]. Dechesne et al. [27] proved that using an energy minimization framework is feasible for smoothing patch maps. We consider our results to form a good basis for incorporating into such label-smoothing approaches.

The delineation results are typically evaluated by comparing automatically delineated stands with manually generated ones [18,24,27] in addition to using statistical consistency metrics. Such comparisons were not possible in our study due to the lack of available ground truth data. Yet, manually delineated forest stand maps have been claimed to be subjective rather than representing a reliable reference [29], and stand delineation in general has been described as a somewhat arbitrary process [24]. Therefore, we do not consider our evaluation approach to be a major limitation to the findings, and we consider the consistency validation to be meaningful.

Nevertheless, a comparison of our results to ground truth data is pending in order to evaluate the accuracy of the delineated forest patches. For delineations that search for patches with homogeneous ecological process properties, as in Experiment 2, it would be conceivable to also incorporate measurements taken at multiple sample points in the field as reference data.

A more fundamental aspect of patch delineation that became apparent in our study is the semantic meaning of forest type delineations in general. The two experiments that were carried out illustrate that the delineation task has multiple solutions rather than a single, unique one. In particular, the delineation of areas with a similar canopy height structure differs remarkably from the one that aims to detect forest patches with similar water cycle conditions (Figure 7 vs. Figure 13; Figure 8 vs. Figure 14, respectively). This finding emphasizes that, firstly, patch delineation is to be seen in regard to a specific task. Secondly, the general concept of patch delineation must also be reconsidered. Patch delineation aims to detect similar canopy structures in order to define distinct class breaks. However, particularly in natural forest areas with a diverse structure or in cases where multiple forest type classes coexist within a certain area [19], canopy properties are continuous in space, with no distinct breaks. Segmentation, however, aims to categorize a continuous canopy property into a reduced number of classes which, in our view, cannot be absolute and permits a certain range of possible solutions.

5.3.2. Definition of Number of Patch Classes

Defining a stopping criterion for segmentation is an important aspect for adequate forest patch delineation and relates to the detection of under- or oversegmentation. In the case of segmenting point clouds to objects, e.g., trees as in Amiri et al. [55], a model of the object under study can be introduced in order to determine the adequacy of the level of the segmentation.

However, for forest patch delineations, such a model cannot be formulated. The evaluation of segmentation has long been a topic in computer vision where, for unsupervised segmentation, within-segment homogeneity should be large, and segments should differ significantly from neighboring regions (e.g., Zhang et al. [56]). Metrics measuring within-segment homogeneity and variations in-between them were applied to forest stand delineations in Mustonen et al. [34] and Tokola et al. [57]. However, those metrics capture the accuracy of the segmentation rather than defining a stopping criterion for the segmentation. Therefore, Wu et al. [29] introduced a minimum size above which forest patches are retained, while segments whose sizes fall below this threshold are merged with neighboring ones until all segments reach the minimum size. Koch et al. [24], on the other hand, defined a number of forest classes into which the forest should be delineated. Such an approach, however, requires information on expected forest classes.

In this study, we intended to avoid pre-defining expected classes, as the aim was to establish a purely data-driven patch delineation approach. That is, we intended to delineate forest patches with similar conditions solely based on the structures detected in the feature space for a set of attributes.

Comparing the consistency of the segmentation with the resulting segment sizes (Table 4 and Table 8 for the Burgenland scene; Table 5 and Table 9 for the Ötscher scene) can hint at the appropriate number of forest patch classes. For the Burgenland scene, the most significant increase in consistency is achieved when cluster merging is stopped at $k=5$ instead of $k=4$ classes (increase of $R^2=0.05$ and $R^2=0.07$ for h95 and $vg\_layerDiff$, respectively), while the increase in consistency is smaller when stopping the merging process of overlapping clusters at $k=8$ instead of $k=5$ classes (increase of $R^2=0.01$ for both attributes). For the Ötscher scene, the largest increase in consistency is also achieved when stopping the merging step at $k=5$ instead of $k=4$ classes (increase of $R^2=0.11$ and $R^2=0.08$ for h95 and $vg\_layerDiff$ attributes), while the consistency increase equals zero when stopping the merging step at $k=7$ instead of $k=5$ classes. Reducing the number of classes from $k=5$ to $k=4$, therefore, marks a crucial point for both scenes while, for larger numbers of patch classes, the consistency comes to a certain plateau.

However, it is up to the user to define what number of classes is appropriate for the intended task. This also has to be considered in light of the resulting class sizes that decrease below a meaningful size to be reasonably retained, if too many classes are retained. Our experiments reveal that, although no explicit minimum class size is defined in the last merging step, merging overlapping patch classes primarily sets in at the smallest clusters, which is an advantage when tackling the size issue.

5.3.3. Sensitivity to Thresholds

The number of resulting forest classes in the framework is not directly determined by a value but results from the similarity criteria in the splitting step (threshold $th\_split$) and in the merging step (threshold $th\_mergeOverlap$). The definition of the intended number of classes thus amounts to the parameterization of these two thresholds.

Sensitivity to $th\_mergeOverlap$

The analyses in our experiments show that the number of final forest classes is mainly driven by the threshold $th\_mergeOverlap$, which controls the overlap that is tolerated between two clusters in one feature. As depicted in Table 10, we found that the resulting number of forest classes has a nonlinear dependence on $th\_mergeOverlap$. The sensitivity of the number of classes to $th\_mergeOverlap$ falls within a certain range, where larger deltas in $th\_mergeOverlap$ are required to change the number of classes. Moreover, this range of decreased sensitivity to $th\_mergeOverlap$ exactly overlaps with the range of the number of classes that we determined to be reasonable when considering the segmentation consistencies along with the class sizes. The coincidence of reduced sensitivity with the range of reasonable numbers of classes is a clear strength of the approach, and we expect that it is due to the overlap criterion $ftr\_mergeOverlap$ having a physical meaning that relates to forest structure.

Sensitivity to $th\_split$

Experiment 1 in the Burgenland scene reveals clear differences in the final delineations for the class size and shape (Figure 17), as well as for the patch class consistencies (Table 12), depending on the applied value for the threshold $th\_split$. Increasing the number of initial splitting steps increases the inter-segment homogeneity of the initial segments. However, it is likely that the forest scene after this processing step is delineated into too many segments which are of too small a size. Merging such small segments in the current implementation of the framework is based on the analysis of the standard deviation of one attribute only. This is in contrast to the splitting step that performs clustering using the entire provided feature set and, therefore, optimizes the homogeneity with regard to all provided features. Furthermore, with too many splitting iterations, oversegmentation results in entities close to superpoints. While these superpoints reveal a very high within-homogeneity, the drawback is that spatially adjacent superpoints possibly display considerable differences in the feature space and are therefore not merged in the subsequent steps. However, one would reasonably disregard such small-scale dissimilarities in order to achieve a certain level of generalization.

The results of the final delineation when employing too severe initial oversegmentation in the splitting step ($th\_split=0.1$) differs remarkably from the result where an intermediate threshold ($th\_split=0.4$) was used. For fewer initial splitting iterations ($th\_split=1.0$), in contrast, the final result is closer to the intermediate one ($th\_split=0.4$).

Considering the consistency of the resulting delineation depending on the stipulated $th\_split$ value, we found a peak in the consistency for a certain threshold value. As discussed for the parameterization of $th\_mergeOverlap$, $th\_split$ can therefore be tuned while considering the consistency metric $R^2$. Again, its parameterization is specific to both the scene and the attribute set used for the delineation.

6. Conclusions

In this study, we established a forest patch delineation framework based on point entities in 3D ALS data. The framework is based on an iterative k-means bipartitioning that uses a set of features that describe the forest canopy structure on multiple scales. The framework allows for the incorporation of any set of possible features describing the forest characteristics. Furthermore, no prior knowledge of the expected patch characteristics or patch classes within a forest scene is required as the framework is entirely data-driven.

Based on experiments using two test sites with different forest types, compositions, and management types and over different topographies, we can attest the robustness of the implemented method and the validity of the patch delineations. The analysis of resulting patches from two different homogeneity definitions reveals that the delineated patches well reflect the spatial distribution of canopy characteristics and that the patches show high within-class homogeneities (from $R^2=0.60$ to $R^2=0.91$).

We see the main utility of our approach for applications in the context of natural managed forests that are characterized by high vertical diversity and alterations through succession processes. Furthermore, we expect ecological studies that rely on forest patches with similar structural canopy properties as basic units for their analyses to benefit from our segmentation framework, as it allows the user to decide which structural information should be incorporated into the delineation process.