ISPRS Journal of Photogrammetry and Remote Sensing 104 (2015) 44–52 Contents lists available at ScienceDirect ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs Effect of slope on treetop detection using a LiDAR Canopy Height Model Anahita Khosravipour a,⇑, Andrew K. Skidmore a, Tiejun Wang a, Martin Isenburg b, Kourosh Khoshelham a,c a Department of Natural Resources, Faculty of Geo-Information Science and Earth Observation, University of Twente, P.O. Box 217, 7500 AE, Enschede, The Netherlands Rapidlasso GmbH, Friedrichshafener Straße 1, 82205 Gilching, Germany c Department of Infrastructure Engineering, University of Melbourne, Victoria 3010, Australia b a r t i c l e i n f o Article history: Received 16 September 2014 Received in revised form 23 February 2015 Accepted 23 February 2015 Keywords: LiDAR Processing Forestry Tree detection Point cloud a b s t r a c t Canopy Height Models (CHMs) or normalized Digital Surface Models (nDSM) derived from LiDAR data have been applied to extract relevant forest inventory information. However, generating a CHM by height normalizing the raw LiDAR points is challenging if trees are located on complex terrain. On steep slopes, the raw elevation values located on either the downhill or the uphill part of a tree crown are heightnormalized with parts of the digital terrain model that may be much lower or higher than the tree stem base, respectively. In treetop detection, a highest crown return located in the downhill part may prove to be a ‘‘false’’ local maximum that is distant from the true treetop. Based on this observation, we theoretically and experimentally quantify the effect of slope on the accuracy of treetop detection. The theoretical model presented a systematic horizontal displacement of treetops that causes tree height to be systematically displaced as a function of terrain slope and tree crown radius. Interestingly, our experimental results showed that the effect of CHM distortion on treetop displacement depends not only on the steepness of the slope but more importantly on the crown shape, which is species-dependent. The influence of the systematic error was significant for Scots pine, which has an irregular crown pattern and weak apical dominance, but not for mountain pine, which has a narrow conical crown with a distinct apex. Based on our findings, we suggest that in order to minimize the negative effect of steep slopes on the CHM, especially in heterogeneous forest with multiple species or species which change their morphological characteristics as they mature, it is best to use raw elevation values (i.e., use the un-normalized DSM) and compute the height after treetop detection. Ó 2015 Published by Elsevier B.V. on behalf of International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). 1. Introduction Information on individual trees is critical for a variety of forest activities and for environmental modeling at the local and regional scales (Lichstein et al., 2010). In the last decade, airborne Light Detection and Ranging (LiDAR) has become a reliable remote sensing technique for estimating individual tree parameters, due to its capability to generate detailed and very precise three-dimensional tree information (Hyyppä et al., 2008; Lim et al., 2003). As an initial and important step in any analysis of LiDAR data on individual trees, treetop detection has attracted much attention and research (Hosoi et al., 2012; Hyyppä et al., 2012; Jing et al., 2012; Kaartinen et al., 2012; Popescu and Wynne, 2004; Vastaranta et al., 2011). Identifying the correct treetop can provide ⇑ Corresponding author. Tel.: +31 630205659. E-mail addresses: a.khosravipour@utwente.nl (A. Khosravipour), a.k.skidmore@utwente.nl (A.K. Skidmore), t.wang@utwente.nl (T. Wang), martin@rapidlasso. com (M. Isenburg), k.khoshelham@unimelb.edu.au (K. Khoshelham). accurate information on crown characteristics and the tree height information, which in turn are useful inputs for growth and volume estimation models (Gebreslasie et al., 2011; Vastaranta et al., 2011; Wulder et al., 2000). A widespread approach is to identify local maxima, which generally correspond to the location and height of individual trees, and then to construct crown segments (Falkowski et al., 2006; Næsset and Økland, 2002; Solberg et al., 2006; Véga and Durrieu, 2011). The local maxima are typically obtained from the height variation of a LiDAR-derived Canopy Height Model (CHM), also known as a normalized Digital Surface Model (nDSM) (Forzieri et al., 2009; Li et al., 2012; Persson et al., 2002; Yu et al., 2011). There are two ways to create a CHM: with rasters or with point clouds (Li et al., 2012; Persson et al., 2002). When working with rasters, the LiDAR ground returns are used to create a raster DTM (Digital Terrain Model), and the highest or first LiDAR returns are used to create a raster DSM (Digital Surface Model). Then the raster DTM is subtracted from the raster DSM to create the final raster CHM (Lim et al., 2003). When working with point clouds, the http://dx.doi.org/10.1016/j.isprsjprs.2015.02.013 0924-2716/Ó 2015 Published by Elsevier B.V. on behalf of International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). A. Khosravipour et al. / ISPRS Journal of Photogrammetry and Remote Sensing 104 (2015) 44–52 classified LiDAR is height-normalized by replacing the raw elevation of each return (i.e. its z coordinate) with its height above the DTM (Khosravipour et al., 2014; Van Leeuwen et al., 2010). Either way, the end result is the absolute canopy height above the bareearth terrain surface. Although the procedure of computing local maxima from a CHM is conceptually simple, the accuracy of its result largely depends on the quality of the acquired LiDAR data, its processing and/or post-processing, and the forest conditions (Kaartinen et al., 2012). For example, the use of a higher density of laser pulse footprints improves the chance of the laser hitting the treetops (Hyyppä et al., 2008; Lefsky et al., 2002), and the use of an efficient local maxima technique enhances treetop identification by reducing commission and omission errors (Chen et al., 2006; Kaartinen et al., 2012; Vauhkonen et al., 2012). A new study suggests that the accuracy of treetop detection can be improved further by removing height irregularities in the CHM (Khosravipour et al., 2014). Moreover, a number of studies indicate that the various forest conditions (e.g., crown sizes, ages, site types, tree species and forest density) can significantly influence intermediate LiDAR derivatives and thereby the performance of tree detection algorithms (Falkowski et al., 2008; Pitkänen et al., 2004; Popescu and Wynne, 2004; Vauhkonen et al., 2012; Yu et al., 2011). Complex forest terrain presents a challenging problem, as it affects the performance of the height normalization step by distorting the CHM, which can reduce the accuracy of extracted tree biophysical parameters (Vega et al., 2014). On steep slopes, the raw elevation values located, for example, on either the downhill or the uphill part of a tree crown are height-normalized with parts of the DTM that may be much lower or higher than the tree stem base, respectively (Breidenbach et al., 2008). Therefore, in the CHM, the downhill part of the crown will ‘‘rise’’ while the uphill part will ‘‘sink’’, causing the entire tree crown to be systematically distorted. In treetop detection, the ‘‘rising’’ branch overhanging lower terrain in the downhill part can turn into a ‘‘false’’ local maximum that is distant from the true treetop. This problem was posed in Isenburg’s keynote speech at Silvilaser 2012 (Isenburg, 2012). He found a CHM that overestimated true tree height by more than double: eucalyptus trees on steep and eroded slopes in the Canary Island of Tenerife were estimated as being 51 m tall whereas their true height was 25 m. Takahashi et al. (2005) and Véga and Durrieu (2011) also reported that one of the sources of tree height overestimation from LiDAR-derived CHM is a horizontal offset error between field and LiDAR treetop detection, particularly on steeper slopes. They concluded the difference may be due to the LiDAR-derived treetop simply being identified as the maximum value of CHM within the crown area on steeper slopes. Heurich et al. (2003) pointed out that this error increases for leaning trees and/or steeper terrain slopes. Breidenbach et al. (2008) reported an increasing underestimation of the CHM-derived height with steeper upward slopes and vice versa for downward slope, which can cause tree height – one of the most important stand characteristics determined in forest inventory – to be misinterpreted, thereby affecting estimates of subsequent biophysical parameters such as biomass, volume and carbon sequestration. The recent study of Vega et al. (2014) suggested using un-normalized elevation values (i.e. using the DSM), and computing the height after a tree crown segmentation step, to avoid the undesirable effect of steep slopes on the CHM. However, until now, the influence of the normalization process on treetop detection and height estimation has neither been studied nor quantified. The aim of this study was to quantify the effects of slope gradient on the accuracy of treetop detection when using a LiDARderived CHM. We first present a simplified theoretical model to illustrate how normalization causes a systematic error in CHMbased treetop detection when an individual tree is located on a 45 slope. We then assess the accuracy of treetop detection by using both the CHM (i.e. the normalized elevations) and the DSM (i.e. un-normalized elevations). Next, we compute the positional difference between the same tree detected in both the CHM and the DSM, in order to investigate the influence of the slope on the horizontal displacement of CHM-detected trees and its effect on subsequent height estimation. 2. Theoretical model The systematic error in CHM-based treetop identification can be quantified by using a conceptual model that is based on field measurement of tree heights. In the field, the original tree height is determined as a vertical distance from tree apex to the upslope root crown (Husch et al., 1982). According to the model (illustrated graphically in Fig. 1), the height of a tree is calculated as the magnitude (length) of a vector h originating at the base of the tree and ending at the treetop. Without loss of generality, we can assume that the tree height is formulated as: h ¼ b þ 2r ð1Þ where b is the crown base height and r is the radius of the hypothetical tree crown. When computing the tree height from the height-normalized model (i.e., CHM) the distance from the highest crown return to its projection on the DTM is used, which introduces a systematic error when the terrain is sloping (Takahashi et al., 2005; Véga and Durrieu, 2011; Vega et al., 2014). We assume a tree on a terrain of constant slope with an idealized spherical crown is always hit at the highest point across the canopy by the laser pulses (i.e. no Fig. 1. Schematic diagram of the geometry involved in the treetop detection based on the effect of slope gradient on a LiDAR-derived CHM. 46 A. Khosravipour et al. / ISPRS Journal of Photogrammetry and Remote Sensing 104 (2015) 44–52 canopy penetration). Although a tree crown is never have perfectly spherical in nature, tree crowns (e.g., coniferous or deciduous) are typically considered to be circular in nadir view (Biging and Gill, 1997; Doruska and Burkhart, 1994). Then we can use the following simplified theoretical model and estimate the CHM-derived tree height as the local maximum of the function hCHM(x) within r 6 x 6 r: hCHM ðxÞ ¼ b þ cðxÞ þ sðxÞ ð2Þ where b is the constant crown base height, and c(x) and s(x) are the contribution of the tree crown and slope in the estimation of CHMderived tree height, respectively. The crown contribution c(x) includes a constant r as the crown radius, and a vertical distance that is a function of the horizontal displacement x from the original treetop (see Fig. 1). This perpendicular distance can be calculated with the Pythagoras theorem that relates the lengths of the three sides of any right triangle as: cðxÞ ¼ r þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r 2  x2 ð3Þ The slope contribution s(x) describes a vertical distance that is also a function of the horizontal displacement x from the original treetop. It is based on the terrain slope that is assumed to be constant below the tree: sðxÞ ¼ mx ð4Þ where m is the constant terrain slope. The resulting function hCHM(x) expresses the tree height that we expect to measure for our idealized circular crown along a line that goes through the true treetop in the direction of the steepest slope for a given radius r and a constant slope m. The xmax is the x that locally maximizes the function hCHM(x) within the crown diameter (i.e. r 6 x 6 r). The xmax estimates the expected systematic horizontal displacement from the true treetop (x = 0). The value hCHM(xmax) estimates the expected CHM-derived tree height. The two extrema xext – one local maximum and one local minimum – of the function hCHM(x) can be determined by finding the zero crossings of its derivative function hCHM0 (x) so that the xmax can be found by solving hCHM0 (x) = 0: hCHM ðxÞ ¼ b þ r þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r 2  x2  mx ð5Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 hCHM ðxÞ ¼ x= r 2  x2  m pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 hCHM ðxext Þ ¼ x= r 2  x2  m ¼ 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xext ¼ mr= m2 þ 1 The horizontal displacement is always directed downhill. Represented in our schematic diagram (Fig. 1), the horizontal displacement will be negative for positive slope (m > 0) but positive for negative slopes (m < 0). horizontal displacement ¼ xmax pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xmax ¼ mr= m2 þ 1 ð6Þ By substituting the x value in the CHM-derived tree height function (Eq. (2)) with the result obtained for xmax, we can estimate the expected CHM-derived tree heights. The expected error difference (vertical displacement) can be calculated by subtracting the true tree height (Eq. (1)) from the expected CHM-derived one (Eq. (2)): vertical displacement ¼ hCHM ðxmax Þ  h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  m2 þ 1  1 vertical displacement ¼ r ð7Þ Though the theoretical model over-simplifies crowns in the real world, it is adequate to demonstrate the systematic error introduced by sloping terrain during the normalization step: the height of LiDAR canopy returns is overestimated downhill of the tree stem and it is underestimated uphill of the tree stem. This systematic error is one of the sources of error that usually leads to tree height being misinterpreted due to steep slope (Breidenbach et al., 2008; Takahashi et al., 2005; Véga and Durrieu, 2011). The formulas in Eqs. (6) and (7) suggest that for a constant slope (m) both the expected horizontal and vertical displacement increase linearly when the crown radius (r) increases (an elliptical crown behaves similarly). However, for a constant radius (r), when slope increases, slope effect on horizontal displacement increases asymptotically whereas slope effect on vertical displacement increases exponentially. As an example, we simulate the theoretical model for different slopes for a constant radius (r) of 3.5 m (Fig. 2). The asymptotic effect on horizontal displacement is clearly visible when the slope approaches 40 degrees. The effect of slope on vertical displacement increases exponentially, especially when the slope reaches 45°, and goes to infinity as the slope approaches 90°. 3. Experimental data 3.1. Study area Bois Noir (44°230 N, 6°450 E) is a forest in the Barcelonnette basin in the southern French Alps (Fig. 3). The northern part of Bois Noir is characterized by irregular rugged topography with slope gradients ranging from 10° to 35° (Thiery et al., 2007) and the southern part is characterized by extensive steep scree slopes of up to 70° (Razak et al., 2011b). The test area was 1.30 km2, with a cover predominantly of mountain Pine (Pinus uncinata) and Scots pine (Pinus sylvestris), at an elevation ranging from 1400 to 2380 m above sea level (Razak et al., 2013). 3.2. Field data Field inventory data were collected during September 2011 and 2012. Some parts of the forest had been affected by a landslide that had caused the tree stems to bend and tilt. We used the landslide boundaries provided by Thiery et al. (2007) and Razak et al. (2011a) to establish 46 plots outside the landslide area. For this research, the measurements collected included tree location, tree crown diameter (in two perpendicular directions) and tree species determination (Table 1). The position of individual trees and the central point of each plot were recorded using the Leica 1200 Differential GPS System and a total station (see Khosravipour et al. (2014) for more detail). The total number of trees sampled was 514. 3.3. LiDAR data and pre-processing The LiDAR data were acquired using a helicopter in July 2009. The small footprint full-waveform LiDAR system (RIEGL VQ-480i) utilized by Helimap has been developed specifically for mapping mountainous forested area (Vallet and Skaloud, 2004). The system was operated at a laser pulse repetition rate of 300-kHz and a scan width of 60° and performed on-line full-waveform analysis to extract up to five discrete returns for each pulse. The survey was flown 250 m above ground level, resulting in a mean footprint diameter of 75 mm on the ground. In order to increase the laser pulse density, the area was covered by seven overlapping flight lines. The mean point density of all return was 160 points/m2. The LiDAR data were stored separately in adjacent, non-overlapping tiles. The LiDAR points were retiled to a tile size of 300 m with a 15 m buffer along all sides of each tile, in order to avoid edge artifacts at the tile boundaries during tile-based processing (Isenburg et al., 2006b). The LiDAR points were classified A. Khosravipour et al. / ISPRS Journal of Photogrammetry and Remote Sensing 104 (2015) 44–52 47 Fig. 2. The relationship between slope, horizontal displacement, and vertical displacement for an idealized spherical crown with a radius of 3.5 m. Fig. 3. Location of the Barcelonnette basin in the map of France (left) and the slope map of the Bois Noir forest (right). into ground and non-ground returns, using automatic progressive TIN (triangular irregular network) densification filtering as developed by Axelsson (2000). Once classified, the ground returns point clouds were interpolated with a TIN by Delaunay triangulating (Isenburg et al., 2006a), which was then rasterized onto a grid with 0.15 m2 spatial resolution, to create the DTM. This LiDAR-derived DTM of the Bois Noir forest was quantitatively and qualitatively assessed, using procedures described by Razak et al. (2011b). The vertical accuracy of the DTM varied between 0.28 and 0.36 m compared to the field data, depending on whether the terrain was open or forested. From this DTM a slope gradient raster was computed with ESRI’s ArcGIS software, which implements the third-order 48 A. Khosravipour et al. / ISPRS Journal of Photogrammetry and Remote Sensing 104 (2015) 44–52 Table 1 Descriptive statistics of the tree crown diameter measurements (m). Minimum Maximum Median Mean Std Dev All (n = 514) Scots pine (n = 263) Mountain pine (n = 251) 0.50 9.10 2.50 2.82 1.47 0.90 9.10 3.45 3.66 1.40 0.50 6.70 1.70 1.94 0.92 finite difference method (Horn, 1981). The slope at each tree stem location used in our analysis was calculated as the mean terrain slope within the reference tree crown. The first returns were used for generating the DSM and the CHM. Such LiDAR-derived surface models often contain so-called ‘‘data pits’’ which occur, for example, when a laser pulse penetrates deep into the canopy before producing its first return, or when multiple flight lines that scan the canopy from different viewpoints are merged (Axelsson, 1999; Ben-Arie et al., 2009; Leckie et al., 2003). The pit-free algorithm, developed by Khosravipour et al. (2014), was used for generating the pit-free CHM. This algorithm consists of two stages: the first stage normalizes the height of the LiDAR returns and generates a standard CHM raster from all first returns and several partial CHM rasters from only those first returns that are above a set of increasing height thresholds (e.g., above 2 m, 5 m, 10 m and 15 m). The second stage composes the standard and all partial CHM rasters into one final CHM by keeping for each pixel the highest value across all CHMs. A variation of the pit-free algorithm was used to generate the pit-free DSM simply by omitting the normalization step and operating on the original elevations. The CHM and DSM were rasterized to the same 0.15 m2 grid spacing as that of the DTM. The pre-processing was implemented by batch-scripting several modules of the LAStools software (rapidlasso GmbH, 2014). 3.4. Individual treetop detection In order to extract individual treetops, a method based on morphological opening and reconstruction was applied to both the CHM and the DSM. The mathematical morphological opening operation (erosion followed by dilation) with an appropriate structuring element (which defines a neighborhood around a given point) is an image-processing technique that can separate different objects (Serra, 1982; Vincent, 1993) and preserve the structural information of each object, based on the structuring element’s shape (Wang et al., 2004). The shape and size of the structuring element are commonly based on the shape and size of the objects of interest. In this study, a flat disk with a diameter of 1.05 m (7 pixel sizes) was experimentally found appropriate (based on the field crown size data) to identify the treetops. The opening operations removed ‘‘foreground’’ objects (i.e., treetops) that were smaller than the selected structuring element in the image. The result is an opened image. Afterward, morphological reconstruction, which is a very efficient method for extracting regional maxima (Khoshelham et al., 2010; Vincent, 1993), was implemented. The opened image was selected as a marker under the original canopy surface as a mask image, in order to retrieve the shape of tree crown boundaries that were smoothed by the opening operation. Subsequently, the reconstructed image was subtracted from the original canopy surface in order to isolate the individual treetops that had been removed by the opening operation as regional maxima. The local maxima of each regional component were extracted from the image. They are the estimated treetop points (x, y, and z). For the CHM, the z coordinate of each treetop point corresponds to the estimated tree height. For the DSM, the z coordinate of each treetop point was height-normalized using its projection onto the DTM. 3.4.1. Accuracy assessment of individual treetop detection The performance of the treetop detection was evaluated by comparing the automatically detected trees from both the CHM and DSM with the trees measured in the field. If one treetop had been detected within a reference crown boundary, the detection was considered correct. If more than one treetop was detected, the closest reference tree was considered as correct, and other trees were then defined as commission errors. However, if no LiDAR-detected treetop was found, this error was considered as an omission error. Subsequently, the distance between the trees detected from the CHM and those detected from the DSM was measured, in order to calculate the horizontal displacement of detected trees in the CHM. The linear regression model was used to find the relationship between the horizontal displacement and the slope of the terrain. Due to the lack of accurate tree height measurements in the field for our super-high density LiDAR points, it was not feasible to validate the tree heights derived from the CHM and the DSM. Thus, only the differences between the height measurements in the CHM and the DSM were used to find the relationship between the vertical displacement and the slope of the terrain. In addition, a visual 3D comparison between the CHM and the DSM provided further insight into the results. 4. Experimental results In order to create a CHM for detecting individual treetop position, the normalization process was applied to the elevation of the LiDAR points. Fig. 4 shows the distorting influence of slope on this LiDAR height normalization on the original morphological structure of the crown of a Scots pine and the crown of a mountain pine. As can be seen, the Scots pine tree, which has a wider, more irregular crown pattern and weaker apical dominance, is affected more than the mountain pine tree, which has a smaller and narrower crown. In order to quantify the effect of the slope, the treetop detection technique was applied to both the CHM and DSM. Fig. 5 shows the CHM raster, examples of trees detected in both the CHM and the DSM, and also omission and commission errors. The Fig. 5 also illustrates an example of the horizontal displacement of a CHM-detected treetop position compared to its location as detected in the corresponding DSM. Table 2 presents the numbers and the percentages of correctly detected trees, omission and commission errors, assessed for Scots pine and mountain pine. The overall tree detection rate of both species was high (‘‘Correct’’ Table 2). Note that the number of false trees results (‘‘Commission’’) was higher for Scots pine than for mountain pine. Out of a total 514 trees, 395 trees were detected in both the CHM and the DSM. The range of positional differences between the same CHM-detected and DSM-detected trees was between 0.15 m (one pixel) and 1.80 m. Because remotely sensed data of high spatial resolution is often spatially auto-correlated (Chen et al., 2012; Hu et al., 2014), we used the range of semi-variogram statistics in order to adjust the minimum distance of the horizontal displacement of treetops. The semi-variograms (using the spherical model) of the both CHM and the DSM images indicated that pixels are highly correlated within a range of 0.35 m. Therefore, treetops detected in both the CHM and DSM within 0.35 m of each other were assumed to belong to the same tree (i.e. to have a positional difference of 0.0 m). The numbers and percentages of the same treetops that were correctly detected from both the CHM and the DSM are reported in Table 3 per tree species and for three different slope ranges (i.e. 0°–15°, 15°–30°, and 30°–45°). On a terrain slope of less than 15° there was no significant horizontal displacement of treetops 49 A. Khosravipour et al. / ISPRS Journal of Photogrammetry and Remote Sensing 104 (2015) 44–52 Fig. 4. The effect of slope on the LiDAR data before normalization (a) and after normalization (b) for Scots pine and mountain pine on a slope gradient of approx. 35°. Fig. 5. An example of identified treetops, omission and commission errors, and the horizontal displacement between a CHM-detected and a DSM-detected treetop position. Table 2 Tree detection results for the CHM and the DSM. Field-measured trees CHM DSM Number of trees Species Correct n (%) Omission n (%) Commission n (%) Correct n (%) Omission n (%) Commission n (%) 263 251 Scots pine Mountain pine 215 (81.7) 207 (82.4) 48 (18.2) 44 (17.5) 30 (11.4) 1 (0.3) 223 (84.8) 206 (82.0) 40 (15.1) 45 (17.9) 25 (9.5) 0 (0.0) 514 Total 422 (82.1) 92 (17.8) 31 (6.0) 429 (83.4) 85 (16.5) 25 (4.8) 50 A. Khosravipour et al. / ISPRS Journal of Photogrammetry and Remote Sensing 104 (2015) 44–52 Table 3 Percentages and numbers of correctly detected trees located on three different slopes. Field-measured trees Number of trees Slope class Linked CHM- and DSM-detected trees Species Total n Without displacement n (%) With displacement n (%) 64 179 20 0–15 15–30 30–45 263 81 148 22 251 0–15 15–30 30–45 Horizontal displacement (m) Vertical displacement (m) Min Max Mean Min Max Mean Scots pine 44 142 15 43 (97.7) 126 (88.7) 8 (53.3) 1 (2.2) 16 (11.2) 7 (46.6) 0.42 0.42 0.80 0.42 1.71 1.80 0.74 1.39 0.10 0.03 0.01 0.10 0.75 1.78 0.16 0.97 Total 201 177 (88.0) 24 (11.9) 0.42 1.80 0.92 0.01 1.78 0.40 Mountain pine 64 114 16 64 (100) 111 (97.3) 16 (100) 0 (0.0) 3 (2.6) 0 (0.0) 0.45 1.27 0.73 0.07 0.14 0.10 Total 194 191 (98.4) 3 (1.5) 0.45 1.27 0.73 0.07 0.14 0.10 between the CHM and the DSM, but the number of affected trees increased with steepness of the terrain. The maximum displacement of CHM-detected treetops from DSM-detected treetops was 1.80 m and occurred on slopes of more than 30°. As expected, the horizontal positional error caused an increasing vertical displacement error (i.e., an overestimation of CHM-derived tree height). The maximum height overestimation was 1.78 m for slope of more than 30°. However, the effect of the systematic error became evident only for Scots pine trees, not for mountain pine trees. A remarkable 46.6% of correctly detected Scots pine trees located on steep terrain (more than 30°) were affected by the systematic error, whereas the conical mountain pines were not affected by this error. A linear regression model using slope as the independent variable explained 54% of the variance associated with the positional displacement of CHM-detected Scots pine trees (Fig. 6). The result of its ANOVA demonstrated that the regression model is significant when predicting the horizontal displacement of CHM-detected trees (F = (1,23) = 23.73, p < 0.05). The field-measured crown radius of Scots pine species was also plotted against the horizontal displacement. Surprisingly, however, there was no correlation between the magnitude of the systematic error and the crown size within Scots pine trees (Fig. 7). Fig. 7. Horizontal displacement of the Scots pine treetops regressed against crown radius. The height differences between the CHM-derived tree height and DSM-derived tree height were plotted against slope gradient. On average, the CHM-derived height was 0.40 m above the DSMderived height, as was expected from our theoretical model. The absolute minimum difference was 0.01 m and absolute maximum was 1.78 m. The regression model using slope terrain explained 48% of the variance associated with the vertical displacement rate (Fig. 8). Its ANOVA result demonstrated the regression model is a significant model when predicting the height differences by slope (F = (1,23) = 20.44, p < 0.05). 5. Discussion Fig. 6. Horizontal displacement of the Scots pine treetops regressed against slope. One of the challenges when generating a CHM through height normalization is that a systematic error appears whenever trees are located on complex, sloping terrain. Based on this observation, we theoretically and experimentally quantified the effect of slope on the accuracy of treetop detection and its effect on the tree height estimation before and after height normalization. The theoretical model presented the systematic error in LiDARderived treetop detection and height estimation as a function of terrain slope surface and tree crown radius. For a constant radius, A. Khosravipour et al. / ISPRS Journal of Photogrammetry and Remote Sensing 104 (2015) 44–52 Fig. 8. Vertical displacement of the Scots pine trees’ height regressed against slope. the horizontal displacement increases asymptotically with the slope and reaches its maximum at a value equal to the radius size; the vertical displacement increases exponentially and goes to infinity as the slope approaches 90°. For a constant slope, both the horizontal (i.e., treetop) and vertical (i.e., height) displacement increase linearly, concomitantly with the crown radius. However, our idealized circular tree crown is different from the crowns that were measured as an average of two longest extending perpendicular directions in the field. This explains why we did not observe any correlation between the magnitude of the systematic error and the field crown size measurement. To better approximate other crown shapes, it would be possible to replace the circular crown contribution to the CHM-derived tree height estimate with various non-circular functions. Our idealized theoretical model is intended to show that terrain slope systematically creates an error for all crown shapes but the actual magnitude of this error will vary with each individual shape. Our experimental results showed that the effect of CHM distortion on treetop displacement depends not only on the steepness of the slope but – more importantly – on the crown shape, which is mainly determined by species. The influence of the systematic error was only evident for the Scots pine trees (which are have a larger average crown size), not for the mountain pine trees (which have a narrow, conical crown with distinct apex). Holmgren and Persson (2004) reported that tree species conical in shape, such as the mountain pine, have a lower 90th LiDAR height percentile than Scots pine. This suggests that the LiDAR returns from the apical portion of mountain pine (i.e. the highest 10%) are distinctly higher than other canopy returns. Hence, no matter how steep the slope, the apical portion always becomes the local maximum in both the CHM and the DSM and is always identified as the correct treetop. Our experimental results also showed that the effect of systematic error on horizontal displacement varies greatly among Scots pine trees. The probable reason is that the crown shape of Scots pine tends to differ with the age of the tree. Scots pine trees are conical when young but become more rounded and irregular as they mature (Holmgren and Persson, 2004; Ross et al., 1986). Local tree competition can also influence crown shape, with the resulting more slender and conical crown (Rouvinen and Kuuluvainen, 1997) being similar to the crown shape of mountain pine. The higher percentage of the commission error also 51 confirmed that the Scots pine crown shape is more irregular than that of mountain pine. Previous studies such as Popescu and Wynne (2004) and Kato et al. (2009) have reported that rounder crowns resulted in individual trees being falsely delineated by local maxima techniques. According to the theoretical model, the height of LiDAR returns calculated from the distance between the highest crown return and its projection on the DTM on steep areas will be overestimated. Our experimental results confirmed that the average vertical displacement error between CHM-derived tree height and DSM-derived tree height was a positive value, i.e., an overestimation of tree height by CHM. The cause was obviously the horizontal downhill displacement of CHM-detected trees on sloping terrain. A few researchers have suggested that the tree height overestimation error might be caused by the differences between field-measured and LiDAR-derived heights on steeper slopes (Breidenbach et al., 2008; Takahashi et al., 2005; Véga and Durrieu, 2011). Our results demonstrate for the first time that – in steep terrain – the height normalization step causes systematic vertical errors, due to the horizontal displacement of treetops. Furthermore, we have shown that these errors depend on the crown shape of the species being measured. On the basis of our results, we recommend using raw elevation values (i.e., the un-normalized DSM) and computing the height after treetop detection and/or tree crown segmentation, in order to minimize the negative effect of steep slopes on the CHM, especially a heterogeneous forest type consisting of multiple species, or species which change their morphological characteristics as they mature. 6. Conclusions We have theoretically and experimentally shown in this paper that the LiDAR height normalization process systematically reduces the accuracy of treetop detection when trees are located on sloping terrain. Our experiments also showed that the effect of slope-distorted CHMs on treetop detection strongly depends on the tree’s crown shape, which is largely determined by its species. It would be interesting to investigate the effect of height normalization on various other non-conical crown shapes (e.g., flat or ellipsoidal) located on terrain with a more complex topography (e.g. with gullies, mounds and other sudden local elevation changes). Acknowledgments The authors would like to express their gratitude to Menno Straatsma for delivering high quality LiDAR data. We also would like to thank Joy Burrough for editing the language of a near-final draft of the paper. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.isprsjprs.2015.02. 013. References Axelsson, P., 1999. Processing of laser scanner data—algorithms and applications. ISPRS J. Photogr. Remote Sens. 54, 138–147. Axelsson, P., 2000. DEM generation from laser scanner data using adaptive TIN models. Int. Arch. Photogr. Remote Sens. 33, 110–117. Ben-Arie, J.R., Hay, G.J., Powers, R.P., Castilla, G., St-Onge, B., 2009. Development of a pit filling algorithm for LiDAR canopy height models. Comput. Geosci. 35, 1940– 1949. Biging, G.S., Gill, S.J., 1997. Stochastic models for conifer tree crown profiles. Forest Sci. 43, 25–34. 52 A. Khosravipour et al. / ISPRS Journal of Photogrammetry and Remote Sensing 104 (2015) 44–52 Breidenbach, J., Koch, B., Kändler, G., Kleusberg, A., 2008. Quantifying the influence of slope, aspect, crown shape and stem density on the estimation of tree height at plot level using lidar and InSAR data. Int. J. Remote Sens. 29, 1511–1536. Chen, Q., Baldocchi, D., Gong, P., Kelly, M., 2006. Isolating individual trees in a savanna woodland using small footprint lidar data. Photogr. Eng. Remote Sens. 72, 923–932. Chen, Y.X., Qin, K., Liu, Y., Gan, S.Z., Zhan, Y., 2012. Feature modelling of high resolution remote sensing images considering spatial autocorrelation. Int. Arch. Photogr. Remote Sens. Spatial Inform. Sci. 1, 467–472. Doruska, P.F., Burkhart, H.E., 1994. Modeling the diameter and locational distribution of branches within the crowns of loblolly pine trees in unthinned plantations. Can. J. For. Res. 24, 2362–2376. Falkowski, M.J., Smith, A.M.S., Hudak, A.T., Gessler, P.E., Vierling, L.A., Crookston, N.L., 2006. Automated estimation of individual conifer tree height and crown diameter via two-dimensional spatial wavelet analysis of lidar data. Can. J. Remote Sens. 32, 153–161. Falkowski, M.J., Smith, A.M.S., Gessler, P.E., Hudak, A.T., Vierling, L.A., Evans, J.S., 2008. The influence of conifer forest canopy cover on the accuracy of two individual tree measurement algorithms using lidar data. Can. J. Remote Sens. 34, S338–S350. Forzieri, G., Guarnieri, L., Vivoni, E.R., Castelli, F., Preti, F., 2009. Multiple attribute decision making for individual tree detection using high-resolution laser scanning. For. Ecol. Manage. 258, 2501–2510. Gebreslasie, M.T., Ahmed, F.B., Van Aardt, J.A.N., Blakeway, F., 2011. Individual tree detection based on variable and fixed window size local maxima filtering applied to IKONOS imagery for even-aged Eucalyptus plantation forests. Int. J. Remote Sens. 32, 4141–4154. Heurich, M., Schneider, T., Kennel, E., 2003. Laser scanning for identification of forest structure in the Bavarian forest national park. In: Proc. Scandlaser Scientific Workshop on Airborne Laser Scanning of Forests, Umeå, pp. 97–106. Holmgren, J., Persson, Å., 2004. Identifying species of individual trees using airborne laser scanner. Remote Sens. Environ. 90, 415–423. Horn, B.K.P., 1981. Hill shading and the reflectance map. Proc. IEEE 69, 14–47. Hosoi, F., Matsugami, H., Watanuki, K., Shimizu, Y., Omasa, K., 2012. Accurate detection of tree apexes in coniferous canopies from airborne scanning light detection and ranging images based on crown-extraction filtering. J. Appl. Remote Sens. 6. Hu, B., Li, J., Jing, L., Judah, A., 2014. Improving the efficiency and accuracy of individual tree crown delineation from high-density LiDAR data. Int. J. Appl. Earth Obs. Geoinf. 26, 145–155. Husch, B., Miller, C.I., Beers, T.W., 1982. Forest Mensuration, third ed. John Wiley & Sons Inc., New York. Hyyppä, J., Hyyppä, H., Leckie, D., Gougeon, F., Yu, X., Maltamo, M., 2008. Review of methods of small-footprint airborne laser scanning for extracting forest inventory data in boreal forests. Int. J. Remote Sens. 29, 1339–1366. Hyyppä, J., Yu, X., Hyyppä, H., Vastaranta, M., Holopainen, M., Kukko, A., Kaartinen, H., Jaakkola, A., Vaaja, M., Koskinen, J., Alho, P., 2012. Advances in forest inventory using airborne laser scanning. Remote Sens. 4, 1190–1207. Isenburg, M., 2012. The story of LAStools. In: Proc. Silvilaser 2012, Vancouver, 19 September, Keynote speech. Isenburg, M., Liu, Y., Shewchuk, J., Snoeyink, J., 2006a. Streaming computation of delaunay triangulations. Proc. ACM Trans. Graph. 25, 1049–1056. Isenburg, M., Liu, Y., Shewchuk, J., Snoeyink, J., Thirion, T., 2006b. Generating raster DEM from mass points via TIN streaming. In: Raubal, M., Miller, H., Frank, A., Goodchild, M. (Eds.), Geographic Information Science. Springer, Berlin Heidelberg, pp. 186–198. Jing, L., Hu, B., Li, J., Noland, T., 2012. Automated delineation of individual tree crowns from LiDAR data by multi-scale analysis and segmentation. Photogr. Eng. Remote Sens. 78, 1275–1284. Kaartinen, H., Hyyppä, J., Yu, X., Vastaranta, M., Hyyppä, H., Kukko, A., Holopainen, M., Heipke, C., Hirschmugl, M., Morsdorf, F., Næsset, E., Pitkänen, J., Popescu, S., Solberg, S., Wolf, B.M., Wu, J.-C., 2012. An International comparison of individual tree detection and extraction using airborne laser scanning. Remote Sens. 4, 950–974. Kato, A., Moskal, L.M., Schiess, P., Swanson, M.E., Calhoun, D., Stuetzle, W., 2009. Capturing tree crown formation through implicit surface reconstruction using airborne lidar data. Remote Sens. Environ. 113, 1148–1162. Khoshelham, K., Nardinocchi, C., Frontoni, E., Mancini, A., Zingaretti, P., 2010. Performance evaluation of automated approaches to building detection in multi-source aerial data. ISPRS J. Photogr. Remote Sens. 65, 123–133. Khosravipour, A., Skidmore, A.K., Isenburg, M., Wang, T., Hussin, Y.A., 2014. Generating pit-free canopy height models from airborne lidar. Photogr. Eng. Remote Sens. 80, 863–872. Leckie, D., Gougeon, F., Hill, D., Quinn, R., Armstrong, L., Shreenan, R., 2003. Combined high-density lidar and multispectral imagery for individual tree crown analysis. Can. J. Remote Sens. 29, 633–649. Lefsky, M.A., Cohen, W.B., Parker, G.G., Harding, D.J., 2002. Lidar remote sensing for ecosystem studies. Bioscience 52, 19–30. Li, W., Guo, Q., Jakubowski, M.K., Kelly, M., 2012. A new method for segmenting individual trees from the lidar point cloud. Photogr. Eng. Remote Sens. 78, 75– 84. Lichstein, J.W., Dushoff, J., Ogle, K., Chen, A., Purves, D.W., Caspersen, J.P., Pacala, S.W., 2010. Unlocking the forest inventory data: relating individual tree performance to unmeasured environmental factors. Ecol. Appl. 20, 684–699. Lim, K., Treitz, P., Wulder, M., St-Onge, B., Flood, M., 2003. LiDAR remote sensing of forest structure. Prog. Phys. Geogr. 27, 88–106. Næsset, E., Økland, T., 2002. Estimating tree height and tree crown properties using airborne scanning laser in a boreal nature reserve. Remote Sens. Environ. 79, 105–115. Persson, Å., Holmgren, J., Soderman, U., 2002. Detecting and measuring individual trees using an airborne laser scanner. Photogr. Eng. Remote Sens. 68, 925–932. Pitkänen, J., Maltamo, M., Hyyppä, J., Yu, X., 2004. Adaptive methods for individual tree detection on airborne laser based canopy height model. Int. Arch. Photogr. Remote Sens. Spatial Inform. Sci. 36, 187–191. Popescu, S.C., Wynne, R.H., 2004. Seeing the trees in the forest: using lidar and multispectral data fusion with local filtering and variable window size for estimating tree height. Photogr. Eng. Remote Sens. 70, 589–604. Rapidlasso GmbH, 2014. LAStools Rapid LiDAR Processing. <http://www. rapidlasso.com> (accessed 30.01.14). Razak, K.A., Buksch, A., Damen, M.C.J., van Westen, C.J., Straatsma, M.W., de Jong, S., 2011a. Characterizing tree growth anomaly induced by landslides using LiDAR. In: Proc. 2nd World Landslide Forum, Rome, pp. 235–241. Razak, K.A., Straatsma, M.W., van Westen, C.J., Malet, J.P., de Jong, S.M., 2011b. Airborne laser scanning of forested landslides characterization: terrain model quality and visualization. Geomorphology 126, 186–200. Razak, K.A., Bucksch, A., Straatsma, M., Van Westen, C.J., Bakar, R.A., De Jong, S.M., 2013. High density airborne lidar estimation of disrupted trees induced by landslides. In: Proc. International Geoscience and Remote Sensing Symposium, Melbourne, pp. 2617–2620. Ross, J., Kellomäki, S., Oker-Blom, P., Ross, V., Vilikainen, L., 1986. Architecture of Scots pine crown: phytometrical characteristics of needles and shoots. Silva Fennica 20, 91–105. Rouvinen, S., Kuuluvainen, T., 1997. Structure and asymmetry of tree crowns in relation to local competition in a natural mature Scots pine forest. Can. J. For. Res. 27, 890–902. Serra, J., 1982. Image Analysis and Mathematical Morphology. Academic Press Inc., London. Solberg, S., Naesset, E., Bollandsas, O.M., 2006. Single tree segmentation using airborne laser scanner data in a structurally heterogeneous spruce forest. Photogr. Eng. Remote Sens. 72, 1369–1378. Takahashi, T., Yamamoto, K., Senda, Y., Tsuzuku, M., 2005. Estimating individual tree heights of sugi (Cryptomeria japonica D. Don) plantations in mountainous areas using small-footprint airborne LiDAR. J. Forest Res. 10, 135–142. Thiery, Y., Malet, J.P., Sterlacchini, S., Puissant, A., Maquaire, O., 2007. Landslide susceptibility assessment by bivariate methods at large scales: application to a complex mountainous environment. Geomorphology 92, 38–59. Vallet, J., Skaloud, J., 2004. Development and experiences with a fully-digital handheld mapping system operated from a helicopter. Int. Arch. Photogr. Remote Sens. Spatial Inform. Sci. 35 (Part B), 791–796. Van Leeuwen, M., Coops, N.C., Wulder, M.A., 2010. Canopy surface reconstruction from a LiDAR point cloud using Hough transform. Remote Sens. Lett. 1, 125–132. Vastaranta, M., Holopainen, M., Yu, X., Hyyppä, J., Mäkinen, A., Rasinmäki, J., Melkas, T., Kaartinen, H., Hyyppä, H., 2011. Effects of individual tree detection error sources on forest management planning calculations. Remote Sens. 3, 1614– 1626. Vauhkonen, J., Ene, L., Gupta, S., Heinzel, J., Holmgren, J., Pitkänen, J., Solberg, S., Wang, Y., Weinacker, H., Hauglin, K.M., Lien, V., Packalén, P., Gobakken, T., Koch, B., Næsset, E., Tokola, T., Maltamo, M., 2012. Comparative testing of single-tree detection algorithms under different types of forest. Forestry 85, 27–40. Véga, C., Durrieu, S., 2011. Multi-level filtering segmentation to measure individual tree parameters based on Lidar data: application to a mountainous forest with heterogeneous stands. Int. J. Appl. Earth Obs. Geoinf. 13, 646–656. Vega, C., Hamrouni, A., El Mokhtari, S., Morel, J., Bock, J., Renaud, J.P., Bouvier, M., Durrieu, S., 2014. PTrees: a point-based approach to forest tree extraction from lidar data. Int. J. Appl. Earth Obs. Geoinf. 33, 98–108. Vincent, L., 1993. Morphological grayscale reconstruction in image analysis: applications and efficient algorithms. IEEE Trans. Image Anal. 2, 176–201. Wang, L., Gong, P., Biging, G.S., 2004. Individual tree-crown delineation and treetop detection high-spatial-resolution aerial imagery. Photogr. Eng. Remote Sens. 70, 351–357. Wulder, M., Niemann, K.O., Goodenough, D.G., 2000. Local maximum filtering for the extraction of tree locations and basal area from high spatial resolution imagery. Remote Sens. Environ. 73, 103–114. Yu, X., Hyyppä, J., Vastaranta, M., Holopainen, M., Viitala, R., 2011. Predicting individual tree attributes from airborne laser point clouds based on the random forests technique. ISPRS J. Photogr. Remote Sens. 66, 28–37.