Impact of Tree-Oriented Growth Order in Marker-Controlled Region Growing for Individual Tree Crown Delineation Using Airborne Laser Scanner (ALS) Data

Abstract

Region growing is frequently applied in automated individual tree crown delineation (ITCD) studies. Researchers have developed various rules for initial seed selection and stop criteria when applying the algorithm. However, research has rarely focused on the impact of tree-oriented growth order. This study implemented a marker-controlled region growing (MCRG) algorithm that considers homogeneity, crown size, and shape using airborne laser scanning (ALS) data, and investigated the impact of three growth orders (i.e., sequential, independent, and simultaneous) on tree crown delineation. The study also investigated the benefit of combining ALS data and orthoimagery in treetop detection at both plot and individual tree levels. The results showed that complementary data from the orthoimagery reduced omission error associated with small trees in the treetop detection procedure and improved treetop detection percentage on a plot level by 2%–5% compared to ALS alone. For tree crown delineation, the growth order applied in the MCRG algorithm influenced accuracy. Simultaneous growth yielded slightly higher accuracy (about 2% improvement for producer’s and user’s accuracy) than sequential growth. Independent growth provided comparable accuracy to simultaneous growth in this study by dealing with overlapping pixels among trees according to crown shape. This study provides several recommendations for applying region growing in future ITCD research.

1. Introduction

Individual tree measurement is an important component of forest inventory that was historically dependent on highly trained personnel in expensive and time-consuming field surveys [1]. Since the application of aerial photographs to forest inventory in the mid-20th century, the efficiency of field work has been complemented by image analysis. However, obtaining individual tree-based information through visual interpretation of imagery is still a cost- and labor-intensive process [2,3]. Thus, researchers have devoted significant effort to developing techniques to obtain individual tree- and stand-level information by automatically delineating individual tree crowns in remotely sensed data, (e.g., [4–6]). The application of very high spatial resolution imagery and small footprint LiDAR (light detection and ranging) data supports these modern approaches to forest inventory [7]. Accurately delineated tree crowns are an essential component of precision forestry to support estimation of crown size and tree diameter [8], crown closure [9], canopy structure [10], height and biomass [11], and improve tree species classification [7] and tree growth evaluation [12].

Improvement in the spatial resolution capabilities of remote sensing instruments has played a vital role in the advancement of new applications [13]. Higher spatial resolution data has stimulated research into semi- and fully-automatic individual tree crown delineation (ITCD) in the forestry, remote sensing, and computer vision communities over the past two decades. Automated ITCD usually requires data with ground sample distance <1 m. One recent trend in the field is the application of active data sources and the integration of high spatial resolution active and passive data. The predominant active remotely sensed data used in ITCD is small-footprint airborne laser scanning (ALS) data from laser range measurements (LiDAR). Three-dimensional ALS data provide a large number of statistical features (e.g., tree height) with high spatial resolution that are difficult to generate from multispectral imagery. Hyyppä et al. [14] used high-pulse-rate laser scanners to detect single trees in boreal forest zones in order to retrieve stand attributes, such as mean height, basal area, and stem volume. Yu et al. [12] used small footprint, high sampling density ALS data to detect harvested trees and estimate tree growth. Forzieri et al. [15] presented a new time- and cost-effective procedure to automatically detect tree position and estimate crown boundaries and plant density using a canopy height model (CHM) with 1 m pixel size derived from ALS data. Researchers also expect benefits by combining information from multispectral imagery and LiDAR. Leckie et al. [5] isolated individual trees using LiDAR and multispectral imagery separately, and then combined delineated crowns from the two datasets to estimate tree height at both individual tree and stand levels. The integration of passive and active data also benefits tree species classification and other individual tree parameter estimation (e.g., [14,16]), and has shifted the focus of some researchers from ITCD algorithm development to extraction of forest parameters. However, the benefit of combining LiDAR and multispectral imagery in ITCD procedure is still not clearly defined [5].

Researchers have developed a variety of semi- and fully-automated algorithms for ITCD. The basic assumption of most ITCD studies is that a treetop is located at a point with maximum radiometric (multispectral imagery) or height (ALS) value and these values decrease toward the crown boundary (e.g., [5,17,18]). Many algorithms, such as marker-controlled watershed segmentation [19] and the tree identification and delineation algorithm (TIDA), developed by Culvenor [20], require treetop detection, e.g., through local maximum filtering, as a step prior to crown delineation. Region growing is another classic image processing method used in tree crown delineation that typically requires treetops as an initial input. Region growing starts from a set of seed points and grows regions until some stop criteria are satisfied [21]. There are two significant questions for region growing: (1) where to start; and (2) when to stop. Region growing, in some studies (e.g., [20,22]), uses a “pouring” approach where stop criteria are defined by boundaries determined by other algorithms, such as a watershed or valley-following method. In such cases, the growth is to label pixels within a region rather than define the extent of the region. This study uses the more general definition of region growing that group’s pixels or subregions into larger regions using stop criteria that are inherently part of the growth procedure. Researchers have used different seeds and defined variable rules to control growing procedures [2,23–25] but few ITCD studies consider the impact of growth order of the regions inside the algorithm. Mehnert and Jackway [26] found that tree-oriented growth order may impact accuracy when using a region growing method for ITCD. They sought to resolve the order dependence issues during pixel processing by enhancing Adams and Bischof’s [27] seeded region growing algorithm. However, the impact of order dependence is still not fully resolved.

The objectives of this study are to investigate: (1) the impact of combining ALS data and orthoimagery on treetop detection accuracy; and (2) the impact of growth order (i.e., sequential, independent, and simultaneous) on the accuracy of tree crown delineation in ALS data using a marker-controlled region growing algorithm that considers homogeneity, crown size, and shape of tree crowns. This work will provide recommendations for applying the region growing algorithm in future ITCD studies.

2. Materials

2.1. Study Area

The study area is located within Heiberg Memorial Forest, which is a 1600 ha property owned by the State University of New York College of Environmental Science and Forestry (SUNY-ESF) in Tully, NY (42.75°N, 76.08°W). The terrain at Heiberg includes rolling hills with elevations ranging from 382 m to 625 m above mean sea level. Heiberg is managed to provide plots that represent the forest ecosystems of the Northeastern United States [28]. Deciduous trees mainly comprise red maple (Acer rubrum), sugar maple (Acer saccharum), red oak (Quercus rubra), and various beech (Fagus) and birch (Betula) species; conifers include red and white pine (Pinus resinosa and Pinus strobus), Norway spruce (Picea abies), hemlock species (Tsuga), northern white cedar (Thuja occidentalis) and tamarack (Larix laricina) [29]. This study was conducted in one 1,000 × 1,000 m block at Heiberg Memorial Forest (Figure 1). Two plots of even-aged Norway spruce were selected from within the block: Plot 1 is 1.4 ha (95 × 150 m) and Plot 2 is 2.4 ha (120 × 200 m).

2.2. Remotely Sensed Data

The active sensor data used in this research was discrete multi-return ALS data acquired on 10 August 2010, using an airborne ALS60 sensor flown by Kucera International Inc. [17]. Kucera used the TerraSolid software suite to process the raw laser data and create a canopy point cloud and a bare earth surface [17]. The characteristics of the ALS data are listed in

Table 1. Characteristics of the airborne laser scanning (ALS) data used in this study.

Parameter	Value
Pulse rate	183.8 kHz
FOV	28 degrees
Altitude	487 m
Flying speed	150 knots
Point density	>7 pts/m2 (average: 12.7 pts/m2)

.

Orthoimagery was downloaded from the New York State GIS Clearinghouse ( http://gis.ny.gov/gateway/mg/2006/cortland/) [30]. The four band (blue, green, red, and near infrared) digital aerial orthoimagery had 0.6 m pixel size and was acquired prior to leaf emergence in April 2006. The seasonality was not a concern as the project focused on delineating coniferous crowns. The process used to ensure registration between the orthoimagery and ALS data is described in Section 3.4.4.

2.3. Reference Data

Reference data for the study came from two sources—field measurement and manual image interpretation—and were used during preprocessing, treetop detection, region growing rule development, and accuracy assessment. The field-based dataset (ref_field) includes 770 trees measured in 2008 [31]. Four characteristics were recorded for each tree: stem location, diameter at breast height (DBH), height, and average crown diameter based on measurements in N-S and E-W directions. The descriptive statistics of tree height and average crown diameter in the ref_field dataset are shown in

Table 2. Descriptive statistics of tree height and average crown size for regression.

Statistics	Minimum	Mean	Median	Maximum	Standard Deviation
Tree height (m)	15.4	27.0	27.3	32.8	2.4
Average crown size (m)	1.4	4.2	4.0	10.3	1.4

. The field data were used to establish relationships between various tree characteristics needed to define window size during preprocessing and treetop detection.

An additional dataset of reference trees in the study area block (see Figure 1) was also generated to support the rules used in the region growing algorithm. The ref_ALS dataset included 205 conifers manually delineated from the CHM image of the study block. A set of circles with 20 m radius spaced on a 50 m grid was applied to control the sampling procedure and improve the generalization of the nonprobability sampling [32]. While it is ideal to select several reference trees from each circle, the reference trees needed to be isolated from other trees as much as possible as the ref_ALS dataset was used to establish a relationship between treetop and tree height at the largest crown extension. Thus, not every circle had samples as some had no conifers or no clear conifers that could be used with high confidence.

A dataset of reference treetops (ref_treetop) and reference crowns (ref_crown) in Plot 1 and Plot 2 were used for accuracy assessment. The crowns were visually delineated from the CHM derived from the ALS data and the orthoimagery. The treetops were initially estimated based on the centroids in the ref_crown dataset and then manually adjusted to match the local maximum height in the CHM. These data included all treetops and crowns in Plot 1 and Plot 2 and provided the best means to assess the quality of the detection and delineation results.

3. Methods

3.1. Overview

The methodology used for ITCD in this study included five main steps: data preprocessing, treetop detection, post processing, tree crown delineation, and accuracy assessment. All steps were programmed and implemented in Matlab R 2011b. A flow chart illustrating these steps is shown in Figure 2.

3.2. Data Preprocessing

The purpose of data preprocessing in the study was to derive raster height surfaces from the ALS data. ITCD studies using ALS data typically utilize a CHM or digital surface model (DSM) and assume treetops are points with local maximum height value. Chen et al. [19] found that using a canopy maximum model (CMM) reduced commission error during treetop detection. The CMM considers maximum height in a local neighborhood to avoid pseudo-treetops caused by branches. As summarized in Figure 2, data preprocessing in this study included the following steps: (1) apply inverse distance weighting to interpolate the digital elevation model (DEM) and DSM from first and last return ALS points, respectively; (2) derive a CHM by subtracting the DEM from the DSM; (3) fill data pits in the CHM using a semi-automated pit-filling algorithm [33]; and (4) calculate a CMM using a variable-sized local neighborhood. The window size for creating the CMM was variable and sized to be smaller than the smallest crown size at a specific height value to prevent “valleys” from being filled [19]. For this study, the window size was defined using the 95% lower prediction limit based on a nonlinear regression between tree height and crown size from the ref_field dataset. Gaussian filtering was applied to smooth the generated CMM since not all non-treetop local maxima can be removed for large tree crowns. Gaussian filtering is a common method for suppressing irrelevant local maxima in treetop detection with flexibility in the standard deviation applied [19,34]. The 2-dimensional isotropic Gaussian function can be calculated using Equation (1), where x is the distance from the origin along the horizontal axis, y is the distance from the origin along the vertical axis, and σ is the standard deviation of the Gaussian distribution. In this study, the convolution kernel shown in Equation (2) was applied as a discrete approximation of a Gaussian function with standard deviation of 1.

$$G(x,y)=\frac{1}{2\pi {\sigma }^2}\times e^{-\frac{x^2+y^2}{2{\sigma }^2}}$$

(1)

$$\frac{1}{273}\left[\begin{array}{ccccc}1& 4& 7& 4& 1\\ 4& 16& 26& 16& 4\\ 7& 26& 41& 26& 7\\ 4& 16& 26& 16& 4\\ 1& 4& 7& 4& 1\end{array}\right]$$

(2)

3.3. Treetop Detection

The focus of treetop detection in this study was to identify treetops in the smoothed CMM using local maximum filtering (LMF) with variable window size. As with the initial smoothing of the CMM, the window size was determined by the 95% lower prediction limit of the nonlinear regression between tree height and crown size [19]. The nonlinear regression between these two variables was computed using Equation (3) within the Statistical Analysis System (SAS) 9.3 [35] software package. The nonlinear regression can be transformed to linear regression via a natural-log transformation. If the regression curve was directly used to define the window size, trees with height smaller than the average value would be omitted. Thus, the 95% lower prediction limit was used to reduce omission error in treetop detection. For a regression model Y = Xβ + ɛ, the lower prediction limit is defined by Equation (4), where Ỹ_X is the lower prediction limit at a given tree height x, s² is the mean squared error, t is the inverse of the Student’s t cumulative distribution function with n − 2 degrees of freedom, and S is the covariance matrix of the coefficient estimates, (X′X)⁻¹·s².

$$\mathit{Crown\; size}=e^{0.075+0.048\times \mathit{TreeHeight}}$$

(3)

$$Y_X=Y-t(1-\alpha ;n-2)\cdot \sqrt{s^2+x\cdot S\cdot x^{\prime }}$$

(4)

The relationship between tree height and crown size is illustrated in Figure 3. The solid line is the nonlinear regression curve that represents the average of crown size prediction, while the dashed line is the 95% lower predict limit for a specific height value. While we could reduce α to further reduce omission error, a smaller window would likely produce larger commission error. Thus, α = 0.05 was applied in this study to balance omission and commission errors. Treetops detected in this procedure were used as markers in the marker-controlled region growing algorithm.

3.4. Post-Processing Using Orthoimagery

Using the smoothed CMM, derived from the ALS data, reduced commission error in tree detection that resulted from pseudo-treetops. However, small trees with lower height are poorly defined in the smoothed CMM and are easily omitted, even using the 95% lower prediction limit as the window size. These small trees may have sufficient spectral information to be detected in multispectral imagery, thus, the available orthoimagery was used to complement the smoothed CMM for enhancing treetop detection. Registration of the two data sources was implemented using a correlation coefficient by searching for the best match between the smoothed CMM and orthoimagery [36]. The correlation coefficient of two images is defined by Equation (5),

$$r=\frac{{\sum }_{i=1}^n\left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}{\sqrt{{\sum }_{i=1}^n{\left(X_i-\overline{X}\right)}^2\sqrt{{\sum }_{i=1}^n{\left(Y_i-\overline{Y}\right)}^2}}}$$

(5)

where X_i is the i^th pixel in the smoothed CMM, X̄ is the mean of the smoothed CMM; and Y_i is the i^th pixel of the orthoimage, and Ȳ is the mean of the smoothed orthoimage.

The green band is commonly used in ITCD studies (e.g., [5,31]) and preliminary testing supported the use of the green band for this study. Thus, the same LMF with variable window size used on the smoothed CMM was also implemented on the green band of the orthoimagery. As crowns delineated from multispectral imagery and a CHM are consistent for most trees [5], the local maxima of height and spectral response are expected to be consistent enough to detect the same trees. Thus, for each treetop location detected in the green band, a window was defined by the crown size predicted using Equation (3) and applied to search for a corresponding treetop in the smoothed CMM. For any potential tree identified in the orthoimage that was not matched in the CMM, local height was considered to reduce commission errors due to inclusion of grass or shrubs. If the height value of a potential tree was greater than those of surrounding pixels in a window corresponding to crown size, it was accepted as an omitted tree. Five height thresholds were tested—50th, 55th, 60th, 65th, 70th, and 75th percentiles—meaning that the potential tree was accepted as an omitted tree if its height was greater than 50%, 55%, 60%, 65%, 70%, or 75% of the pixel values in the corresponding crown-sized window.

The majority of ITCD algorithms rely on processing of a single image band. Researchers have applied various bands and derivatives, e.g., principal component analysis, to represent the spectral information in multispectral imagery for ITCD (e.g., [18]). The first principal component (PC1) was tested in the post-processing procedure of this study instead of the green band. PC1 derived from the four band orthoimage carried 92.68% of the information for this block and would appear to be sufficient for ITCD study. However, while PC1 and the green band provided similar results, the green band tended to yield slightly better performance in this study. The benefit of directly applying multispectral bands or integrating multi-sensor data may be explored for ITCD studies on both production and input levels in future. However, registration issues are a significant challenge when multi-sensor data are used.

3.5. Crown Delineation

3.5.1. Marker-Controlled Region Growing

Treetops detected in the previous step were used as initial seeds for marker-controlled region growing. Six conditions, based on homogeneity, crown area, and crown shape, defined the stop criteria for the region growing algorithm. Unlike many region growing algorithms that stop growth when one criterion is satisfied, the six criteria in this study interacted with different focus, including controlling shape of region (Criteria 1 and 2) and growing neighbors (Criteria 3–6). The six criteria considered were: (1) the rectangularity (R) of each region; (2) the ratio of length and width for the region (R_lw); (3) if a pixel was a neighbor of a seed or new starting pixel (NSP), where NSPs are the new neighbors that are added to region in this loop and can grow in the next loop; (4) the standard deviation of height for each region; (5) the region area; and (6) the difference between tree height and height at the largest crown extension. The largest crown extension is the widest part of the tree, that is, point A in Figure 4. The logic of applying these criteria within the growing procedure is shown in Figure 5. Each criterion and the thresholds used are described further below.

Criterion 1: Rectangularity

The shape measure called rectangularity was calculated to control shape of each growing region. Rectangularity (R) is the ratio of the region area to the area of the minimum bounding rectangle and varies between 0 and 1; R is 0.785 (π/4) for a circle. In this region growing algorithm, R for each region was required to be between 0.5 and 1 to avoid regions becoming non-circular.

Criterion 2: Ratio of length and width

It is insufficient to consider only rectangularity since an elongated region may have a relatively rectangular shape, thus a large R. Therefore, the ratio of length and width (R_lw) was also used to control the shape of the growing region with a maximum R_lw threshold of 2 used in this study (R_lw = 1 for a circle). The selection of R_lw threshold was based on trials for including more pixels to grow a crown and avoid a rectangular shape. Criterion 1 and 2 worked together to check region shape at the beginning of every loop (Figure 5). If either Criterion 1 or 2 failed, the growing procedure stopped for this seed.

Criterion 3: Neighborhood

The third stop criterion checked whether a pixel is a neighbor of a seed or NSP. The neighborhood used in this study was the second order rook nearest neighbor with 12 neighbors. The circularity of this neighborhood represented the typical tree crown shape better than alternatives, such as a queen neighborhood. The growing potential of the 12 neighbors was prioritized first based on distance and then on height difference between the seed or NSP and its neighbor, with pixels added to the region sequentially. Therefore, the nearest neighbor with height closest to a seed has the highest probability to grow. Prior to adding a pixel to a region, the remaining stop criteria are checked (Figure 5).

Criterion 4: Threshold of variability

Homogeneity of height within the crown was controlled using the standard deviation of height within each growing region. Researchers have found that the characteristic parameters of image variograms are closely related to forest canopy structures [37,38]. Therefore, the threshold of variability was determined by the variogram of the smoothed CMM used to delineate tree crowns. Figure 6 shows the empirical variogram fitting of the smoothed CMM for two plots using exponential and spherical models, with the exponential model more closely fitting the data. The two variograms did not exhibit the nugget effect (nugget = 0) associated with random variance [37], which indicated that there was no random effect for the height of pixels in these two plots. The overall height variance of the smoothed CMM was 22 m² for Plot 1 and 16 m² for Plot 2. With 0.5 m pixel size, there was no spatial autocorrelation for Plot 1 and Plot 2 in the CMM after 8 m (16 pixels) and 6.5 m (13 pixels), respectively. When a new pixel was added to a tree, the standard deviation of height for the region was calculated. The height standard deviation threshold (thres_std) applied was the square root of variance derived from the variogram based on the maximum distance of pixels in this region. If the standard deviation of height after adding a potential neighbor was greater than thres_std, this neighbor was not added to the region (Figure 5).

Criterion 5: Threshold of crown area

Crown area is a significant factor for tree crown delineation that was used to control overestimation. The area of each growing region should be smaller than the crown area corresponding to a specific tree height. The threshold for crown area (thres_area) was calculated based on a circle with crown size (i.e., diameter) predicted by Equation (3) in the treetop detection procedure. This criterion was checked when a new neighbor is added to the current region for controlling the growing area of each tree crown (Figure 5).

Criterion 6: Threshold of height difference

Theoretically, the difference between tree height and height of any pixel within the crown should be no greater than the difference between the tree height and the height at the largest crown extension (thres_diff). To avoid underestimation, a new pixel with height more than thres_diff was added to region as long as other criteria are satisfied, but this pixel would not have any ability to grow in the next loop. If no further NSP are found among the new neighbors, region growing for this seed stopped (Figure 5).

The thres_diff value was estimated using a voxel-based LiDAR method [39]. A pseudo-waveform showing the frequency for different height values was developed for voxels within each tree (Figure 7a). The point with the highest frequency on the pseudo-waveform was selected to represent the largest crown extension. Pseudo-waveforms were determined for each tree crown in the ref_ALS dataset (described in Section 3.3.3) and height at the largest crown extension (Hgt_crown) was linearly regressed by tree height (Hgt_tree) (Equation (6)) to enable prediction of the largest crown extension for any treetop. The resultant fitting plot with R² of 0.758 is shown in Figure 7b. The height difference threshold (thres_diff) was the difference between these two heights and varied according to tree height.

$${\text{Hgt}}_{\text{crown}}=0.799\times {\text{Hgt}}_{\text{tree}}+2.623$$

(6)

3.5.2. Growth Order in Marker-Controlled Region Growing

Growth order in this study refers to the order that the region growing algorithm completes the growth procedure for all seeds. In this study, three different growth orders were implemented: (1) grow sequentially (G_seq); (2) grow independently (G_ind); and (3) grow simultaneously (G_sim). The common assumption for these three growth procedures, and for most ITCD algorithms, is that pixels can only belong to one tree in the final delineation map. G_seq grows one seed until its stop criteria are satisfied and then grows another seed, completing the growth procedure sequentially (i.e., seed by seed). G_ind grows each seed independently, which results in overlapping pixels among trees after the growth of all seeds. Pixels in any overlapping area are assigned to a single region using a shape factor called circularity (C), based on the assumption of tree crown circularity used in this study. The circularity measure represents how different a region is from a circle and is defined by Equation (7),

$$C=\frac{4\pi \cdot A}{P^2}$$

(7)

where P is perimeter and A is the area of region, and C = 1 for a circle. As this method is only used to assign overlap area after all trees have completed growing, we consider it to be “pseudo-competition”. Figure 8 shows an example where Tree₁ and Tree₂ have an overlap area of O_a (Figure 8a). Pixels in the overlap (O_a) are assigned to a tree based on the circularity (C₁ and C₂) computed. In this case, C₂ is closer to 1 than C₁, thus the overlap area is assigned to Tree₂ (Figure 8b).

The G_sim method is modeled after the idea of an automaton, which is a mechanism to process input information from the surroundings and update the state according to transition rules [40]. Many simple automata have two states representing the presence or absence of a given phenomenon [41]. In this study, each treetop is an automaton that has two states: grow and stop. This automata-based method grows all seeds simultaneously but only grows a portion of the crown at a time, thus, resulting in different growth time depending on crown size. A single growth cycle combines growth loops of NSP for all trees. For a specific cycle, each tree grows NSP, finds NSP for the next cycle, and then waits until all trees have completed the cycle. The stop criteria for G_sim are the same as the G_seq and G_ind methods and spatial movement of crown boundaries takes place through the diffusion of information through neighborhoods [42]. In G_sim, transition rules account for both the state of every automaton and its corresponding neighborhood. A flow diagram of G_sim for region growing is shown in Figure 9.

3.6. Accuracy Assessment

3.6.1. Accuracy Assessment of Treetop Detection

Accuracy assessment is a significant component of ITCD studies and typically has one of two main objectives: (1) determine the accuracy of tree detection, including the number and location of trees (point accuracy); and (2) determine the quality of tree crown delineation, i.e., the boundary of tree crowns (polygon accuracy). This study considered two levels of accuracy assessment of treetop detection and tree crown delineation: plot and individual tree level.

Plot level accuracy (PLA) provides a means to assess the overall performance of a treetop detection procedure while avoiding potential registration issues between reference and delineated trees. PLA is widely used in ITCD research for reporting the quality of tree counts (e.g., [18,31,43]). It is often reported as a detection percentage (DP) by using the ratio of the number of detected trees (N_d) to the number of reference trees (N_r) (Equation (8)). This non-site specific measurement provides an aggregated accuracy for a plot and avoids the need to match tree locations.

$$\mathit{DP}=\frac{N_d}{N_r}\times 100\%$$

(8)

Individual tree level assessment of treetop detection requires evaluation of the accuracy of the location of detected trees. This is a site-specific accuracy in which the detected and reference treetops are compared on a location-by-location basis [44]. As the two points are unlikely to exactly match, detected and reference treetops that match within a given search distance are defined as correctly identified trees. In this study, the 1:1 hit defined by Hirschmugl et al. [24] was applied to assess the accuracy of tree location using a 1 m buffer. A 1:1 hit means that one detected treetop is located within a 1 m buffer of one reference treetop. The small buffer size minimizes the detection of multiple treetops within the buffer of a reference treetop. Producer’s (PA) and user’s accuracy (UA) of treetop detection based on the count of 1:1 hits are defined using Equations (9) and (10), respectively,

$$\mathit{PA}=\frac{N_{1:1}}{N_r}\times 100\%$$

(9)

$$\mathit{UA}=\frac{N_{1:1}}{N_d}\times 100\%$$

(10)

where N_1:1 is the number of 1:1 hits, N_r is the number of reference treetops, and N_d is the number of detected treetops.

3.6.2. Accuracy Assessment of Tree Crown Delineation

The plot level accuracy used to assess the tree crown delineation accuracy in this study was relative error of crown area (RE_CA), which considers the difference between total area of the delineated crowns and that of the reference crowns. It is calculated using Equation (11),

$$\mathit{RE}\_\mathit{CA}=\frac{A_d-A_r}{A_r}\times 100\%$$

(11)

where A_d is the total area of the delineated crowns and A_r is the total area of the reference crowns. A positive RE_CA shows that the ITCD algorithm overestimated total crown area; while a negative value shows that the ITCD algorithm underestimated total crown area.

Researchers have developed various criteria to assess the accuracy of tree crown delineation on an individual tree level. This study used producer’s and user’s accuracy based on 1:1 matched crowns [7,28]. A 1:1 crown match is achieved when the overlap of the delineated and reference crowns is greater than 50% of the area of both of these crowns (Figure 10a). Like treetop assessment, PA and UA for tree crown delineation are defined using Equations (9) and (10), respectively, where N_1:1 is the number of 1:1 matched crowns, N_r is the number of reference tree crowns, and N_d is the number of delineated tree crowns.

Though PA and UA from 1:1 matched crowns provide information about individual tree level accuracy, they do not characterize incorrect delineations. A detailed case analysis of individual tree accuracy was conducted considering both location of treetop and area of tree crown. Nine cases were defined as follows:

(a)

1:1 match—reference crown only overlaps one delineated crown and the overlap area is greater than 50% of both delineated and reference crowns (Figure 10a);

(b)

match but under-grow—reference crown only overlaps one delineated crown and the overlap area is greater than 50% of the delineated crown area but less than 50% of the reference crown area (Figure 10b);

(c)

match but over-grow—reference crown only overlaps one delineated crown and the overlap area is greater than 50% of the reference crown area but less than 50% of the delineated crown area (Figure 10c);

(d)

mis-located match—reference crown only overlaps one delineated crown but overlap area is between 0 and 50% of both delineated and reference crown area (Figure 10d);

(e)

split/multiple match—reference tree is split into multiple delineated trees and at least two delineated trees have overlap area greater than 50% of the delineated crowns (Figure 10e);

(f)

merge—delineated tree merges multiple reference trees and at least two reference trees have overlap area greater than 50% of the reference crowns (Figure 10f);

(g)

multi-intersected—reference tree is intersected by multiple delineated trees and at least one delineated tree has overlap area between 0 and 50% of delineated trees (Figure 10g);

(h)

commission error—delineated tree does not intersect any reference trees (Figure 10h);

(i)

omission error—reference tree does not intersect any delineated trees (Figure 10i).

This analysis simultaneously considers delineated and reference crown perspectives and provides detailed information in terms of merits and drawbacks of the ITCD algorithm.

4. Results

4.1. Impact of Combining ALS Data and Orthoimagery on Treetop Detection

One of the goals of this study was to determine the value in using orthoimagery to complement the ALS data in treetop identification. Trees omitted in the ALS processing were identified using the orthoimagery when they had a local maximum in the green band and were higher than a specific percentile of height within a corresponding crown-size window on the smoothed CMM.

Table 3. Accuracy assessment of treetop detection using ALS and combination of ALS and green band.

Plot	Accuracy Metrics	ALS	ALS & Green Band
			Percentile
			50th	55th	60th	65th	70th	75th
1	DP	80.2%	99.8%	96.4%	93.2%	88.3%	85.5%	82.8%
	N1:1	341	378	378	376	366	363	352
	Nd	377	469	453	438	415	402	389
	Nr	470	470	470	470	470	470	470
	UA	90.5%	80.6%	83.4%	85.8%	88.2%	90.3%	90.5%
	PA	72.6%	80.4%	80.4%	80.0%	77.9%	77.2%	74.9%
2	DP	86.7%	103.1%	97.6%	93.9%	90.2%	88.6%	88.1%
	N1:1	703	723	720	718	715	715	713
	Nd	744	885	837	806	774	760	756
	Nr	858	858	858	858	858	858	858
	UA	94.5%	81.7%	86.0%	89.1%	92.4%	94.1%	94.3%
	PA	81.9%	84.3%	83.9%	83.7%	83.3%	83.3%	83.1%

summarizes the accuracy assessment of treetop detection for Plot 1 and Plot 2 using ALS alone compared to the combination of ALS and the orthoimagery green band. Six percentiles of height in post-processing, 50th, 55th, 60th, 65th, 70th and, 75th, were applied to the detected treetops. For the plot level accuracy, detection percentage decreased (e.g., 99.8% to 82.8% for Plot 1) as percentile increased as DP is a non-site specific accuracy metric and depends only on the number of detected trees. The smaller the percentile used, the more omitted trees identified, and thus the higher DP obtained. In the extreme case, when the 50th percentile was used in Plot 2, the number of detected treetops (885) was greater than the number of reference treetops (858); the 103.1% DP failed to represent the performance of the treetop detection procedure. However, DP still suggests there is value in incorporating the supplemental data as DP was always higher when the two datasets (ALS and green band) were used (80.2% and 86.7% for Plot 1 and Plot 2, respectively), no matter which percentile was considered.

For individual-tree level accuracy, PA and UA of treetop detection using two datasets changed in opposite directions as percentile applied decreased: PAs slightly increased (about 1% to 5%) but UAs decreased (about 10% to 13%). This was because the number of 1:1 hits increased as more treetops were detected, hence PA—the ratio of number of 1:1 hits to the number of reference treetops—increased. On the other hand, as UA was the ratio of the number of 1:1 hits to the number of delineated treetops, it had the potential to decline if the increase in 1:1 hits did not have a corresponding increase in terms of detected treetops. This indicated that not all the detected trees were actual trees, which resulted in commission error increasing with the decrease of percentile. Thus, the 70th percentile was selected as a compromise to balance commission and omission errors for further analysis. This meant that a seed identified in the orthoimagery was included as an omitted tree if its height was larger than 70% of the pixels in its corresponding crown size window. The combination of ALS and green band inputs increased both DP and PA by about 5% and kept a consistent UA. An example of treetop detection using ALS and orthoimagery is illustrated in Figure 11. This figure highlights where a small tree located between two large trees is not initially detected using the CMM, but has sufficient spectral information to be detected using the orthoimagery.

4.2. The Impact of Growth Order on the Accuracy of Tree Crown Delineation

As the combination of ALS and orthoimagery based on the 70th percentile provided the highest detection accuracy, these treetops were used to implement tree crown delineation.

Table 4. Accuracy assessment of tree crown delineation using three region growing methods and the treetops from ALS and green band.

Plot	Accuracy Metrics/Cases Analysis	G_seq	G_ind	G_sim
1	RE_CA	3.8%	1.6%	1.3%
	PA	60.6%	61.5%	62.6%
	UA	70.9%	71.9%	73.1%
	1:1 match	285	289	294
	Match but under-grow	14	11	13
	Match but over-grow	78	80	75
	Mis-located match	13	10	8
	Split/multiple match	0	0	0
	Merge	29	28	29
	Multi-intersected	14	14	13
	Commission error	10	10	10
	Omission error	37	38	38
	No. of Reference Crowns	470
2	RE_CA	−9.0%	−10.2%	−11.1%
	PA	81.1%	82.5%	82.1%
	UA	91.6%	93.2%	92.6%
	1:1 match	696	708	704
	Match but under-grow	13	9	11
	Match but over-grow	20	15	17
	Mis-located match	3	4	4
	Split/multiple match	5	5	5
	Merge	28	30	31
	Multi-intersected	13	7	5
	Commission error	12	12	12
	Omission error	80	80	81
	No. of Reference Crowns	858

summarizes the quantitative assessment of tree crown delineation and shows that region growing using the three different growth methods influenced the delineation results. For plot level accuracy, G_sim tended to reduce total crown area compared to G_ind and G_seq, thus, slightly decreasing the overestimation seen in G_seq (i.e., RE_CA decreased from 3.8% to 1.3%) for Plot 1, but slightly increased the underestimation of G_seq (i.e., decreased RE_CA from −9.0% to −11.1%) for Plot 2. At the individual tree level, G_sim had the best results for plot 1, improving PA from 60.6% (G_seq) to 62.6%, and UA from 70.9% (G_seq) to 73.1%. G_ind had the best results in plot 2, with PA increasing from 81.1% (G_seq) to 82.5% and UA from 91.6% (G_seq) to 93.2%. The accuracy of Plot 2 was generally higher than that of Plot 1 as the trees in Plot 2 had more conical shape and thus were more easily detected and delineated.

Although the three growth methods provided similar results for the two plots, some differences were observed using the nine cases described in Figure 10. For Plot 1, G_sim had the largest number of 1:1 matched trees (294) and G_seq had the smallest number (285); for Plot 2, G_ind had the largest number of 1:1 matched trees (708), while G_sim also had a larger count (704) than G_seq (696). G_seq had slightly weaker results because by completing the growth procedure of every seed sequentially, trees that grow early in the procedure can limit the growth potential of other seeds. The increased growth opportunity of all seeds in G_ind and G_sim resulted in a reduction in the number of “mis-located match” cases (from 13 to 10 or 8) for Plot 1 and a reduction in the number of “multi-intersected” (from 13 to 7 or 5) for Plot 2.

Figure 12 shows an example of tree crown delineation using the marker-controlled region growing with the three growth methods. With sequential growth (G_seq Figure 12a), the delineated crown for tree 1 included part of what should have been tree 2. In some extreme cases, a tree that grows early (e.g., tree 3) can cover the seed of another tree (tree 4) to such a degree that there is no growth potential at all, which resulted in mislocated match and multi-intersected cases for tree 4 (Figure 12a). The independent (Figure 12b) and simultaneous (Figure 12c) approaches avoid this phenomenon, thus, reducing the number of “mis-located match” cases for Plot 1 and the number of “multi-intersected” crowns for Plot 2. For marker-controlled region growing, the main errors in the two plots were “match but over-grow”, “merge”, and “omission error” from omitted treetops.

5. Discussion

The marker-controlled region growing implemented in this study is a two-step ITCD method that requires preliminary treetop detection similar to the approach used in the marker-controlled watershed algorithm applied by Chen et al. [19]. The largest difference between the region growing and watershed approaches is that the marker-controlled watershed algorithm cannot handle gaps between crowns and, thus, a background mask is usually needed if there are large gaps between trees. Even in dense forest, if small gaps exist, the watershed method tends to overestimate the area and perimeter of tree crowns [17], while the marker-controlled region growing approach can handle such gaps without the need for a background mask. This is consistent with the findings of Li et al. [45] that region growing was a better option than other approaches (e.g., valley-following, template matching) in dense forest. The selection of appropriate markers (i.e., treetops) is significant for any marker-controlled ITCD method, and is directly related to commission and omission error. Use of the smoothed CMM in treetop detection avoids commission errors resulting from large branches [19]. Chen et al. [19] found that LMF with variable window size using the 95% lower prediction limit successfully reduced omission error in treetop detection. This study further sought to reduce omission errors by post-processing using orthoimagery as a supplemental dataset. However, the benefit of using multiple datasets (ALS and orthoimagery) greatly depends on successful registration. This study conducted area-based registration based on correlation between the images [36], which is suitable for a forested environment without salient objects [46]. The benefit of using orthoimagery also depended on the percentile of height within a specific crown-size window selected to identify omitted trees. While lower percentiles were associated with detecting more trees, a compromise is needed to balance commission and omission error.

Region growing is a basic segmentation technique that starts from specific seeds and grows homogeneous regions until some stop criteria are satisfied [21]. As with many other ITCD approaches, the marker-controlled region growing algorithm presented in this paper assumes that treetops have the greatest local height in ALS data and these values decrease towards the crown boundary. The algorithm also assumes that trees have a circular shape in two-dimensional imagery, which is suitable for coniferous forests [28]. The algorithm considers homogeneity, shape and region size using six criteria that work in coordination to stop growth of a region. The method of adding new neighbors follows a circular growth trajectory because the neighbors are added first based on distance and then on height difference. Thus, this algorithm delineates homogeneous circle-like crowns according to the height of the seeds (i.e., tree height). Of the criteria used to control region growth, R and R_lw work together to incorporate appropriate pixels in the growth procedure and efficiently control the shape of each crown. The circularity used in the independent growth method is not applied generally to limit growth because the circularity of the second order rook nearest neighbor is too small (<0.5) to start the algorithm. However, circularity was considered a good shape factor to assign overlapping pixels after the growth of each tree. Threshold selection for all of the criteria applied in this algorithm was a significant challenge. This study used linear regression to predict thresholds of crown area (thres_area) and difference between tree height and height at the largest crown extension (thres_diff) to maximize the objectivity of threshold selection and flexibility of the algorithm. As the proposed method incorporates many heuristics based on theoretical advantages, the only way to test the utility of the method is to use real data, especially under varying forest canopy conditions. Kaartinen et al. [47] compared four individual tree detection and extraction techniques using ALS data and found that minimum curvature-based tree detection was a good solution for suppressed trees. Future studies could compare the MCRG method with other algorithms using the same dataset, testing site, and accuracy metrics.

This study demonstrated that growth order can influence the delineation results using marker-controlled region growing. The sequential growth procedure is particularly impacted by seed order. Whatever the order of seeds considered, trees that grow early limit the growth potential of other seeds. Thus, sequential growth yielded the lowest accuracy in this study. The automata-based simultaneous growth method provided slightly higher accuracy than sequential growth because it introduced a time-dimension into the growth procedure and mimics tree growth from an ecological point of view in that trees grow in cycles. However, the automata-based method failed to mimic interaction among trees (i.e., tree competition) because tree growth simply stops if one tree touches other trees. Independent growth also provided comparable or slightly better results than simultaneous growth. However, by using circularity to handle overlapping pixels after the independent growth of all seeds, the independent growth method uses a pseudo-competition related to shape that does not represent natural tree competition. It is likely that this approach may be less suitable when plots have intensive competition because trees with higher circularity will gain all overlapping pixels and limit the growth of neighboring trees. Therefore, while the simultaneous and independent methods had similar accuracy, we believe the simultaneous growth method would provide a more flexible solution with broader applicability. Future studies to develop a competition mechanism that better coincides with natural tree competition may be incorporated into the region growing algorithm to better adjust crown boundaries.

Reference crowns used in accuracy assessment were dominated by manual delineation of crowns from remotely sensed data for most ITCD studies due to the difficulties of obtaining exact tree crowns in fields (e.g., [7,48]). One of the challenges in comparing individual tree crown delineation algorithms across studies is that while many approaches have been suggested, there is still no standardized accuracy assessment procedure. For instance, Leckie et al. [5] defined delineation polygons as isols and reference polygons as refs for assessment using six different cases, a technique that has been used in subsequent studies (e.g., [7]). Jing et al. [48] identified five accuracy categories for assessing tree crown delineation. The 1:1 match case in this study used the same definition of “match” (i.e., overlapping area >50% of both reference and delineated crown) applied by Jing et al. [48]. In order to obtain more details about error, this study split the “nearly matched” scenario used by Jing et al. [48] into “match but under-grow” and “match but over-grow” cases. These were classified into more detailed cases according to overlapping area for better presenting errors resulting from the marker-controlled region growing algorithm. These metrics provide information for adjusting thresholds in the algorithms or for examining threshold impact in future studies.

Ke and Quackenbush [28] also separately considered reference and delineated crown perspectives and reported counts of different ratios of delineated and reference crowns. For simplicity, the delineated: reference crown ratios of 2:1, 3:1, and ≥4:1 from the reference crowns perspective in Ke and Quackenbush [28] were included in the “split/multiple match” case; while 1:2, 1:3, and 1: (≥4) from the delineated crowns perspective were counted in the “merge” case in this study. Therefore, the accuracy assessment in this study combines both delineated and reference perspectives simultaneously for efficiency and convenience. As one delineated crown may be involved in multiple cases (e.g., “mis-located match” and “multi-intersected”), the analysis should be represented by counts of different cases and cannot be converted to a percentage. Producer’s and user’s accuracy related to 1:1 matched trees can be used as percentage accuracy metrics. As 1:1 matched trees are the ultimate goal for ITCD studies, it is better to know not only how many matched trees the producer identifies (producer’s accuracy), but how many trees identified are actually matched (user’s accuracy).

6. Conclusions

This study implemented a two-step marker-controlled region growing algorithm for individual tree crown delineation considering homogeneity, crown size, and shape. The study investigated the benefit of combining ALS data and orthoimagery on treetop detection as well as the impact of growth order (i.e., sequential, independent, and simultaneous) on tree crown delineation. Results showed that the complementary data from the orthoimagery reduced omission error associated with small trees in the treetop detection procedure; however, a compromise in height threshold selection was necessary to balance commission and omission error. Augmentation of reference counts can directly increase PA but may not improve UA if the additional trees are challenging to correctly detect.

This study found that the growth order used in region growing can influence the crown delineation results. The automata-based simultaneous growth yielded the best results with higher accuracy than sequential growth. This is because trees that grow early can limit the growth potential of other seeds with sequential growth. Independent growth also provided a comparable accuracy to simultaneous growth by dealing with overlapping pixels using a pseudo-competition that assigned overlap according to circularity.

As with other ITCD algorithms, marker-controlled region growing has advantages and limitations. The first advantage of this method is avoiding large overestimation of crown area because crown area criteria are checked when a pixel is added to a region, rather than after adding all potential neighbors. Shape criteria are used to control the region structure after adding all new neighbors. The second advantage of this method is processing efficiency and flexibility. It took about five minutes to complete ITCD for the 1.4 ha Plot 1, and seven minutes for the 2.4 ha Plot 2, using Matlab R 2011b on a computer with a 2.3GHz Intel Core i7-3610 QM CPU. The algorithm is also convenient and flexible to use after thresholds are determined by regression or other methods. A further benefit of this region-growing approach is that it does not require preprocessing to eliminate background or non-forest area so is suitable for sites with gaps among trees.

The biggest limitation of the marker-controlled region growing is the selection of treetop markers since the success of treetop detection directly influences the crown delineation results. The complementary orthoimagery provides a reasonable method of finding small trees that would otherwise be omitted using the ALS dataset. The algorithm also requires reliable thresholds in the growing procedure and may have some cost on thresholds selection, e.g., the cost of obtaining field data to establish a relationship between tree height and crown size. The thresholds in this study were selected for Norway spruce. The algorithm will be improved in the future for delineating tree crowns in mixed forest according to different characteristics of various species.

The region growing method used in this study minimally attempts to mimic canopy growth by simply stopping when a growing crown touches neighboring trees. Future ITCD studies should consider tree competition procedures that better conform to the natural process for adjusting overlapping tree crown boundaries. Integration of image processing techniques and ecological processes may become a trend in future ITCD and forest inventory studies.