Workshop ISPRS on Laser Scanning 2007 and SilviLaser 2007, Espoo, September 12-14, 2007, Finland
HIERARCHICAL WATERSHED SEGMENTATION OF CANOPY HEIGHT MODEL FOR MULTI-SCALE FOREST INVENTORY
Kaiguang Zhao *, Sorin Popescu
Spatial Sciences Lab., Dept. of Ecosystem Science and Management, Texas A&M University, College Station, TX 77840, USA - (zhaokg, s-popescu)@tamu.edu
KEY WORDS: LiDAR, CHM, Watershed segmentation, Scale, Crown delineation, Forest inventory
ABSTRACT: Canopy Height Model (CHM) is a standard LiDAR-derived product for deriving relevant forest inventory information, among which individual tree identification is a crucial task. The watershed algorithm from markers is the typical procedure applied to CHMs for delineation of crowns. However, for low-quality CHMs or under certain canopy conditions, segmentation at individual tree level is not practical, e.g., due to grouped trees in dense forests. In this study, we investigated the feasibility of a hierarchical watershed transform (HWT) algorithm to segment CHMs at both individual tree levels and scales above that. As compared to the results by the variable-window filtering for individual trees, HWT allows more flexibilities in removing nontreetop maxima by referring to the dynamic attributes of the potential treetops (i.e., local maxima). It is also found that the choice of filters for smoothing CHM has significant influences on the detection of treetops. Beyond individual tree level, the segmentation by HWT was compared with a commercial package eCognition, and both give similar segmentation results, though with minor differences. Due to the lack of fieldmeasured trees matched with LiDAR-detected ones, no quantitative evaluation of accuracy is provided in this study. Nevertheless, the results of this study reveal that HWT is a viable procedure that could be applied for multilevel segmentation of CHM.
1. INTRODUCTION Reliable mapping of forest resources is a crucial task in many scientific and practical settings, e.g., regional estimate of biomass or fuel models as the input of fire behaviour modelling. For spatially-explicit forest inventory, LiDAR (Light Detection And Ranging) has become a well-established technique in terms of its capability of direct measurements on canopy structures (Hyypp, Inkinen, 1999; Nsset, 2002; Maltamo et al., 2004). Extensive research efforts have been focused on the use of airborne laser scanners for deriving forest information by employing various approaches at relevant analysis unit, i.e., grid, stand or individual tree level. A canopy height model (CHM), which represents the difference between the top canopy surface and the underlying ground topography, becomes a standard LiDAR product that can be effectively derived from LiDAR raw data through appropriate filtering of LiDAR point clouds for the separation of ground hits and canopy hits. In practice, CHMs are available in raster formats and can be considered as 2D images where individual tree crowns are often visually noticeable. To automatically delineate tree crowns or detect individual trees from the CHM, a variety of algorithms or procedures have been devised or explored across various forest conditions, which include but are not limited to image segmentation, local maxima filtering, and template matching (Persson et al., 2002; Popescu et al., 2002; Koch et al, 2006; Chen et al., 2006; Falkowski et al., 2006). Furthermore, with individual trees identified, tree height and crown width also can be measured directly from CHM, and other tree dimension parameters such as stem diameters can be predicted from LiDAR-derived metrics by regression models (Pyysalo and Hyypp, 2002; Popescu et al. 2003). These algorithms for tree identification often make the assumptions that treetops correspond to local maxima present in the CHM, or that tree * Correspondence author
crown assumes certain types of geometry that could be evaluated against the templates used. Among the segmentation approaches, watershed transform is the most popular technique in segmenting a CHM because it is intuitively straightforward to treat each concave tree crown in the inverted CHM as a catchment basin. However, cautions should be exercised as to how to appropriately select local maxima as candidates of treetops. For example, within a single crown, there may be multiple local maxima that result primarily from the real irregularity of crowns or partly from random errors in the procedures of creating the CHM; therefore over-segmentation is usually observed in such situations. As remedies, common strategies are to pre-process the CHM using a smoothing filter, or to merge over-segmented regions ad hoc; but too strong a filter could possibly smear out small trees; as a result, smoothing filters with adaptive parameters are often desired to alleviate such situations. Koch et al. (2006) used a pouringalgorithm, which is similar in spirit to watershed segmentation, to delineate tree crowns over a temperate deciduous and mixed forest, where the authors applied a Gaussian filtering with a varying parameter adaptive to height classes, devised several schemes to deal with spurious regions, and also employed a Ray algorithm to trace the actual crown edge within each segment. On the other hand, in the local filtering for treetops, a smaller window tends to have larger commission errors while a larger window often leads to more omission errors. An improved version of local maxima filtering is demonstrated in Popescu et al. (2003) by using a circular window and adapting its size locally relative to height by referring to a pre-defined height-crown equation; this variable window filtering proves successful considering the fact that higher trees generally have larger crowns. In Chen et al. (2006), a further refinement of variable window size filtering is realized by taking into account the variability in the prediction of crown width from tree height. Falkowski et al. (2006) performed wavelet-transform on CHM
436
IAPRS Volume XXXVI, Part 3 / W52, 2007
image using 2D Mexican Hat wavelet and identified local maxima in the resulting convolved image as potential treetops; their method is in essence a template-matching procedure. Despite all these successes, difficulties or problems are often reported in previous studies on tree-level crown segmentation. These are often witnessed as the relative large error of missing or misidentifying trees, particularly over forests with high canopy closure and density. Ideally, in open forests where no overlapping or suppressed trees exist, all the aforementioned algorithms supposedly are able to find all trees if the shape of tree crowns is also well-defined. Generally, the difficulties in correct delineation may be attributed to the incapability of CHM to capture real canopy surfaces that may be caused by low quality of LiDAR data, improper resolution for rasterization of CHM, or ineffectiveness of procedures for generating CHM (e.g., methods of classifying raw returns, or interpolation algorithms); in most cases, the complexity of canopy surface itself prevents high accuracy in identifying trees. Although sophisticated algorithms tend to consider more aspects either in the phase of pre-processing or post-processing in attempt to reduce commission and omission errors, experiential evidences suggest that in some cases it is extremely difficult or even impossible to delineate individual tree crowns in CHMs even with one of the most sophisticated image processing tools such as the human vision system (Bortolot, 2005). To this end, some researchers, alternatively, went beyond individual tree levels to examine the usefulness of CHM. For example, Bortolot (2006) investigated the use of CHM for tree clusters that correspond to a group of tree crowns. van Aardt et al. (2006) used the eCognition package (Definiens Imaging GmbH, Munich, Germany) to segment CHM at stand levels and then assessed forest volume and biomass on a persegment basis. In certain senses, the use of CHM at scales greater than individual tree levels circumvents the difficulties in crown segmentation and can accommodate the purposes of operational inventory at scales appropriate for forest management. The objective of this study is to investigate segmentation of CHM for forest inventory at multiple scales by using a hierarchical watershed transform algorithm. The watershed algorithm used in this work is a marker-controlled morphological algorithm that has also been utilized for isolating individual trees in previous researches, e.g., Chen et al. (2006), among others. Its hierarchy characteristics result from the use of dynamics as criteria to select markers that are then used for growing basins at the corresponding scale. Specifically, first we rely on the dynamics attributes of local maxima for the selection of potential treetops in the delineation of tree crowns, and next we apply the hierarchical watershed transform (HWT) for CHM segmentation at both tree levels and scales greater than individual trees. The results were compared to those by the established variable window filtering at individual tree levels and those by eCognition at levels above individual trees, respectively.
Forest, many of which with a natural pine stand structure, and upland and bottomland hardwoods. Much of the southern U.S. is covered by forest types similar to the ones included in our study area. 2.2 Airborne Laser Data and Canopy Height Model Laser scanner data were acquired with a Leica-Geosystems ALS40 during the leaf-off season in March 2004 by M7 Visual Intelligence Inc. of Houston, Texas. The LiDAR system was operated to record two returns per pulse, i.e., the first and the last, with a reported accuracy of 20-30 cm and 15 cm for horizontal and vertical positioning, respectively, and was configured to scan +/-10 degrees from nadir. On average, the dataset has a point density of 2.6 hits per m2.
P P
A Digital Surface Model (DSM) was created by first selecting the LiDAR point of maximum height within each 0.5m x 0.5m cell that contains at least one laser hits, and next interpolating the selected laser points into a regular grid with a spatial resolution of 0.5m by triangulation. A Digital Elevation Model (DEM) was derived using a proprietary package and was made available by the data vendor. Consequently, the CHM was generated through the pixelwise subtraction of DEM from DSM. 3. METHODS 3.1 Hierarchical Watershed Transform (HWT) The idea of watershed transform (WS) is typically illustrated with respect to immersion simulation. In classical WS, holes are punched at local minima (to be more precise, regional minima) while a improved algorithm known as watershed from markers (WSM) punches the holes at the prescribed markers (Soille, 2003). Denote the WSM as follows,
R =WSM ( I , M )
where I is the input image, i.e., the inverted CHM, M is the set of markers, and R is the set of segmented regions. Notice that the cardinalities of R and M (the number of elements in each set) are the same and there exists a one-to-one mapping between the two sets. Apparently, if all local minima are considered as makers, the WS and the WSM produce the same results. On the other hand, another algorithm, the Hierarchical watershed transform (HWT), is a multi-scale watershed approach that depends on the dynamic of local minima to create a set of nested partitions (Dougherty and Lotufo, 2003). The dynamic H d of a minimum is defined as the height one has to climb up from the minimum in order to reach another minimum of lower value, as illustrated in Figure 1 (left) for the minimum point m which has a neighbouring minimum m of lower height.. In fact, H d is the height extinction value of the corresponding valley in the h-minima operator; and it has two counterparts, i.e., area-dynamic Ad and volume-dynamic
Vd , which can be defined similarly. For example, the volume-
2. MATERIALS 2.1 Study Area A forested area of approximately 47 sq km, located in eastern Texas, USA, is chosen for this study. The airborne laser coverage consists of pine plantations in various developmental stages, old growth pine stands in the Sam Houston National
dynamic Vd of a minimum is the volume of water that has to be filled to reach another minimum of lower height (Figure 1c). An HWT at a given scale s is the WSM using only local minima with dynamics greater than or equal to s as markers. This can be expressed as,
R = WSM ( I , M s ) = HWT ( I , s)
437
ISPRS Workshop on Laser Scanning 2007 and SilviLaser 2007, Espoo, September 12-14, 2007, Finland
with M s = {m RMIN ( I ); DYN (m) s} where RMIN () is the operator to obtain local minima from the input image I (e.g., the inverted CHM) and DYN () is the operator to calculate the dynamic (or area- and volumedynamics) of a local minimum. When s increases, fewer minima are selected as markers, and hence a coarser segmentation is obtained. Of particular note is that s has no direct correspondence to the scale commonly used for the spatial extents. The implementation of HWT used in this study is based on minimum-cost path algorithm as described in Lotufo and Falcao (2000), and it requires the discrete value of CHM. Thus, we digitized the float height values into integers using a 0.01m quantification interval. If s = 1, i.e., the finest scale, the HWT will produce the same segmentation as classical watershed transform (WS) because all the minima are selected.
area (volume)-dynamic criteria conflict, the latter takes priority as a conservative strategy to reduce over-segmentation. As with (1), the thresholds can also vary adaptively; for example, the threshold for area-dynamic can be the lower limit of predication interval based on a height vs. crown-area equation. The justification for the ratio criterion lies in that the dynamic indicates the depth of crown valley, thus, the deeper the valley is relative to the height of its minimum, the more likely it is to be a crown. However, if the valley is too narrow as indicated by a small value of area-dynamic, it is less possible to be a crown; this argument justifies the area-dynamic criteria. (3) Besides within-crown spurious local minima, there may be nontreetop minima dangling along crown-ground borders that are caused by protruding branches, or in some cases, there exists non-tree features that produce isolated minima. It could be helpful to use H d Ad / Vd (the ratio of dynamic times areadynamic to volume dynamic) as a initial criterion to identify these local minima; for example, if the ratio is near 1, the minimum more likely belongs to this category of nontreetop minima. In this study, the ratio threshold is set to be 0.95. Other more detailed rules could be devised to remove spurious minima or keep treetops. In these rules, adaptive schemes should be preferred if prior knowledge is available, and it also will be advantageous to take into account all the four attributes attached to each minimum. For instance, when using variablewindow filtering, in order not to miss too many treetops, its preferable to have a window size that is a little smaller (i.e., the lower limit of prediction interval of crown width) so as to incorporate the variability of crown width given a tree height; but this leads to a high risk of commission errors. A remedy to alleviate this situation is to refer to Ad or Vd as further guidance. As to the aforementioned rules, of particular note is that for certain local minima, two or more criteria may lead to conflicting judgments; whichever should take precedence is dependent on the degree of belief as to how the assumptions of each criterion approximate the real situations. On the other hand, as another common strategy to reduce local height variations, we pre-processed the CHM by smoothing procedures before applying watershed segmentation. In addition to Gaussian filtering, we also used the wavelet-based filter to de-noise CHM. The use of wavelet for image analysis characterizes the adaptive basis functions for capturing local signal features as well as a multi-scale representation of the image (Matlab Online Help, Mathworks Inc. USA). the Despite availability of automatic de-noising wavelet algorithms with minimal prior input, in this study we adjusted the threshold parameters in wavelet filtering through a trial-and-error approach, as described later in this section. 3.3 Segmentation of CHM beyond individual tree levels Segmentation beyond tree levels is an alternative to analyze CHM when the algorithms of individual tree crown cannot be appropriately applied. As in our case, trees in the CHM over certain forested areas are barely identifiable. To extend the HWT to deal with such cases, there are multiple options for procedures of selecting markers. Each of three types of dynamics, or their combinations could be used as criteria to choose markers from the minima for multi-scale segmentation. For example, in this study, we attempted to choose as markers those minima whose values of H d Ad are larger than a specified threshold; and the threshold plays a role like a scale
m m
m m
m
(a)
(b)
(c)
m
Figure 1. Illustration of the concept of dynamic for a local minimum m whose neighbouring minimum of lower height is
m
; for simplicity, a 1-D signal is used instead of 2D CHM surface. (a) the dynamic of m as indicated by the arrow, (b) the areadynamic of m as indicated by the dashed line, (c) and the volume-dynamic of m as indicated by the hatched area. 3.2 Segmentation of CHM at Individual Tree Level In the segmentation of CHM by watershed from markers, the number of delineated crowns is equal to the number of markers used. Therefore, careful selection of markers as treetops must be performed. It is also impractical to select all the local maxima as treetops due to over-segmentation. Previous researches employed different strategies to perform the subset selection of local maxima (Popescu et al., 2002, Popescu and Wynne, 2004; Chen et al., 2006; Koch et al., 2006). In this work, at least four attributes, i.e., the CHM height and three dynamic values, are tagged to each local minimum. Based on these four attributes directly or other indicators derived from them, a series of decision rules could be devised to help guide the selection of treetop minima in the inverted CHM as demonstrated by a few examples in the following: (1) If the height of a minimum is lower than a threshold, i.e., 2.5 m for this study, it is labelled as nontreetop; however, more complicated schemes could adapt the height threshold locally. The intuition of this rule is that when a minimums height is too low, it is less likely to be a treetop. (2) Given a local minimum with a height larger than a threshold, i.e., 15 m in this work, if the ratio of its dynamic to height is greater than a prescribed value (0.5), it will be identified as treetop; and if its area dynamic or volume-dynamic is below a threshold, it is deemed as a nontreetop. Otherwise, the status of the minimum remains undetermined. In case that the ratio and
438
IAPRS Volume XXXVI, Part 3 / W52, 2007
parameter. Alternatively, volume-dynamic could be directly used as criteria to select markers. As with segmentation of individual trees, smoothing filters can be first applied to CHM for coarser segmentation. Furthermore, with multi-resolution decomposition of CHM by wavelet, we were able to perform segmentation on the wavelet-filtered coarse-level image. In this study, we randomly selected 5 sample subsets of CHM over our study area, each with a size of 256m x 256 m, and applied segmentation to each subset at individual tree levels and above with the aforementioned procedures where we used symlet basis in the wavelet smoothing and decomposition, due to its near symmetry property and its resemblance to crown shape. The 5 selected subsets of CHM represent various growth stages, and all have relatively high canopy closures (e.g., unthinned pine plantations). For the Gaussian filtering, we set = 2 as argued in Chen et al. (2006), and used a window size of 1.5m. In the wavelet-based filtering, we first performed a 2level decomposition of the CHM and then chose leveldependent thresholds for smoothing: at the first level, the threshold was selected as the 90% percentile of magnitudes of detailed coefficients, and at the second level the 70% percentile was used; the thresholds were determined empirically.
left for the smoothed CHM by wavelet filtering. It seems very difficult to recognize individual trees over parts of the CHM.
For the CHM in Figure 2, there are totally 14081 local maxima in the original CHM while the Gaussian filtered CHM only has 3230 maxima as compared to 5305 in the wavelet-filtered CHM. Out of these 14801 local maxima, the variable window filtering (VWF) as proposed in Popescu et al. (2002) identified 2660 of them as treetops. With the criteria using the dynamics properties, 2867 were selected as treetops from the original CHM; with the same dynamic criteria, 1325 local maxima were identified as treetops from the Gaussian filtered CHM, and 2263 from wavelet-filtered CHM. This suggested that the wavelet filter used in this study tends to preserve local features, thus resulting in more local maxima in comparison to Gaussian filtering, as also shown in Figure 2. In all the five selected subsets of CHM, we have a limited number of field-sampled trees. But, unfortunately, for most of these trees, we failed to match them with LiDAR trees detected with the above algorithms. Also, we found it is not an easy endeavour to manually delineate trees out of CHM based on visual interpretation as shown in the close-up view of Figure 2. Therefore, no attempt is made in this study to report quantitatively the accuracy of tree identification due to the unavailability of reference data; and only comparisons between the methods were reported in terms of numbers of detected trees and mean tree height for all the 5 subsets as listed in Table 1. The numbers of detected trees are significantly different among methods (p < 0.005, ANOVA), but the differences in mean tree height are not statistically significant (p= 0.76, ANOVA). Both the two smoothing procedures significantly reduce the tree numbers (p < 0.001, paired-t tests) and the Gaussian filtering produces the least number of trees in all cases. Tree Number
VWF WSOrg WSGau WSWav
4. RESULTS AND DISCUSSION A typical scenario of forested area of the study site was shown in Figure 2 where a portion of the smoothed CHMs respectively by Gaussian and wavelet filters is also displayed, as compared to the original CHM. When evaluated visually it became clear that in our case the Gaussian filtering has stronger smoothing effects than wavelet filter. For example, Gaussian filter can effectively fill the holes within crowns while a certain number of relatively large holes, though reduced, are still preserved in the wavelet-filtered CHM.
Mean Tree height(m)
VWF WSOrg WSGau WSWav
1 2 3 4 5
2660 4468 3583 1416 2572
2867 3910 3419 2355 2720
1325 1587 1504 1302 1229
2263 2702 2597 2221 2017
15.5 12.9 15.2 20.6 14.7
16.0 13.1 15.3 21.8 17.2
15.9 12.2 14.6 22.1 18.8
16.9 13.4 15.5 24.0 20.9
Table 1. Comparison of tree number and mean tree height between different methods where VWF stands for variablewindow filtering in Popescu et al. (2002); and WS-org, WS-gau and WS-wav for watershed segmentation using the dynamicbased criteria applied respectively on original CHM, Gaussianfiltered CHM and Wavelet-filtered CHM.
Figure 2. One selected subset of CHM for this study (above) together with a close-up of the area indicated by the red rectangles (below) where the left is the original CHM, the middle for the smoothed CHM by Gaussian filtering, and the
In addition, qualitative evaluation is given over part of the CHM in Figure 2. It can be seen that no one method is superior to others according to visual examination as demonstrated in Figure 3. However, the smoothing, especially the Gaussian filtering, does help remove some, though not all, spurious local maxima. Also, the smoothing may produce inconsistent results over different parts of the CHM; for example, in Figure 3d there are more trees identified around the centre and fewer trees around the left corner as compared to Figure 3b. Overall, the result for the Gaussian-filtered CHM seems to offer a more satisfactory segmentation on this particular area than other
439
ISPRS Workshop on Laser Scanning 2007 and SilviLaser 2007, Espoo, September 12-14, 2007, Finland
methods, although no optimal selection of filtering parameters and window size was performed.
(a)
(b)
(a)
(b)
(c) (d) Figure 4: Comparison of segmentation results at stand levels: (a) eCognition, (b) HWT on original CHM, (c) HWT on Gaussianfiltered CHM and (d) HWT on the wavelet-filtered CHM 5. CONCLUSION (c) (d) Hierarchical watershed segmentation of CHM is obtained by examining the dynamics properties of local maxima. The use of these dynamic attributes provides extra information as well as more flexibilities in devising rules to determine if a local maximum is treetop or not for individual tree detection. In this study, no sophisticated rules were explored; instead we simply used thresholds for the removal of nontreetop maxima. Further studies could investigate other possible criteria in determination of treetop maxima. Our results also suggested that smoothing plays an important role in suppressing spurious local maxima in CHM, and the Gaussian filter tends to produce stronger smoothing effects than wavelet-based procedures for dense forests of our study area; but neither of the two filters is consistently superior to the other. When it is difficult or infeasible to detect individual from a CHM, HWT is a practical choice to segment CHM at stand level or above. The segm...