Automated extraction and validation of Stone Pine (Pinus pinea L.) trees from UAV-based digital surface models

ABSTRACT Stone Pine (Pinus pinea L.) is currently the pine species with the highest commercial value with edible seeds. In this respect, this study introduces a new methodology for extracting Stone Pine trees from Digital Surface Models (DSMs) generated through an Unmanned Aerial Vehicle (UAV) mission. We developed a novel enhanced probability map of local maxima that facilitates the computation of the orientation symmetry by means of new probabilistic local minima information. Four test sites are used to evaluate our automated framework within one of the most important Stone Pine forest areas in Antalya, Turkey. A Hand-held Mobile Laser Scanner (HMLS) was utilized to collect the reference point cloud dataset. Our findings confirm that the proposed methodology, which uses a single DSM as an input, secures overall pixel-based and object-based F1-scores of 88.3% and 97.7%, respectively. The overall median Euclidean distance revealed between the automatically extracted stem locations and the manually extracted ones is computed to be 36 cm (less than 4 pixels), demonstrating the effectiveness and robustness of the proposed methodology. Finally, the comparison with the state-of-the-art reveals that the outcomes of the proposed methodology outperform the results of six previous studies in this context.


Introduction
Sustainable tree management focuses on maximizing the use of tree resources globally, and various tree species provide socio-economic advantages.Stone Pine (Pinus pinea L.) is a specific coniferous evergreen type of tree wide-spreading, especially in the Mediterranean and Aegean regions (Mutke et al. 2012;Gonçalves and Pommerening 2012;Correia et al. 2018).Among the edible seed forests in the Mediterranean region, Stone Pine is currently the pine species with the highest economic value, with edible seeds called "pine nuts".Given the high economic value of these edible seeds, the Stone Pine stands out among other coniferous forest tree species used for timber harvesting (Schneider, Calama, and Martin-Ducup 2020).This tree is very useful as an industrial product and also provides socio-economic benefits (Buyuksari et al. 2010;Calama et al. 2016;Eker and Laz 2018;Akyol and Örücü 2019).It also has potential as a crop in agroforestry systems in Mediterranean climates, where it has a high value and contributes more to local/national economies.
An issue of strategic importance is quickly and accurately determining physical tree characteristics (e.g.stem location, crown diameter, height, and canopy volume) (Michałowska and Rapiński 2021).
The major focus of early automatic techniques for extracting trees from remotely sensed data was the use of radiometric information in high-resolution 2D images (e.g.Pouliot et al. 2002;Wang, Gong, and Biging 2004;Erikson and Olofsson 2005;Pu and Landry 2012).Despite the limited information available in 2D images, methods solely devoted to the processing of images still represent a major field of tree extraction research (e.g.Jing et al. 2012;Hung, Bryson, and Sukkarieh 2012;Malek et al. 2014;Ozdarici-Ok 2015;Leckie, Walsworth, and Gougeon 2016;Koc-San et al. 2018;Donmez et al. 2021), and receive increased attention due to the recent developments in machine learning, especially in deep learning (e.g.Liu, Wang, and Wang 2019;Wu et al. 2020;Xi et al. 2021).Another large study topic includes experiments whose framework directly incorporates height information from Light Detection and Ranging (LiDAR) data (e.g.Leckie et al. 2003;Hirschmugl et al. 2007;Tooke et al. 2009;Zhen, Quackenbush, and Zhang 2016;Dong et al. 2020).LiDAR technology has been widely used in forestry research for many years (e.g.Michałowska, and Rapiński 2021;Liang et al. 2016;Shao et al. 2020;Campbell 2021).Research promoting either DSM or normalized DSM (nDSM) or a canopy-height model (CHM) is another comprehensive field involving height information.
Unmanned aerial vehicles (UAVs) can now generate highly accurate, detailed, and high-quality data.Because UAVs are usually cost-effective to operate, their usage in forestry is growing rapidly (Chong et al. 2017;Dainelli et al. 2021aDainelli et al. , 2021b)).This paper presents a novel methodology for the automated extraction of individual Stone Pine trees from UAV-based DSMs (Figure 1).First, we present an enhanced probability map of local maxima that uniquely handles the newly proposed local minima information.Next, this new probability map helps us to improve the orientation symmetry evidence used for the detection of regions that are likely to correspond to the Stone Pine trees.Finally, the crown boundaries and stems of individual Stone Pine trees are extracted using the Chan-Vese active contour model (Chan and Vese 2001) that independently processes each detected stem region.
Our main goal in this research is to propose a novel methodology for the automated extraction of Stone Pine trees.Within the framework, as a new contribution, we develop a new enhanced probability map of local maxima that improves the calculation of the orientation symmetry used to detect the tree stem locations.Our presented results not only confirm that the proposed methodology utilizing a single DSM as input is promising and effective but also provide better results than the previous approaches tested in this context.

Related studies
So far, several researchers have particularly studied the automatic extraction of the Stone Pine trees.In the study conducted by Mei and Durrieu (2004), both a LiDAR-based elevation model and multispectral images were jointly used to delineate tree crowns of Pinus pinea in the southern part of France.Pirotti (2010) evaluated a kernel-based approach to extract Pinus pinea tree positions, and the corresponding height using correlation with CHM as a similarity measure in the north-eastern part of Italy.Awad, Jomaa, and Arab (2014)

Study area and preparation of the dataset
Our area of interest is in the Turkish province of Antalya, near the town of Kumkoy in the Aksu district (Figure 2).The whole study area is approximately 10 km 2 of the protected pine forest centered at 30°57ʹ05ʹʹ E and 36°51ʹ59ʹʹ N. The study area is one of the most significant Stone Pine forest locations due to its rich and loamy soils and favorable climatic conditions (mean annual temperature of 29.1°C and mean annual precipitation of 1070 mm).The topography is relatively flat, with a mean height of 3 m above the mean sea level.Although Stone Pine is the most common type of tree in the region, the study area also has a variety of other species such as Red Pine, Bay tree, Aspen tree, and Wild Pear, albeit in small numbers.
An RGB camera (Sony A6000) with a focal length of 30.25 mm and having 24 MP (6000 × 4000 pixels) on board the BRAMOR ppX (GNSS PPK -Post Processing Kinematic) UAV platform was utilized to acquire the images.The UAV was used to undertake several overlapping scans, and the overlap ratios of the UAV images during image acquisition were configured to be high, with 80% forward and 60% side laps.The flight height of the mission was set to approximately 400 m above the ground; thus, it provided images with a ground sampling distance of approximately 5 cm.
The UAV mission began at 11:21 am local time on 16 June 2020 and lasted exactly 40 min.A total of 1492 images were collected during the mission.During imaging, the wind speed was around 18 km per hour, and stratus clouds produced a partial gray layer of cloud cover throughout the study area, causing periods of light precipitation or drizzle (Figure 3).Besides, no GCPs were deployed on the ground before the mission, as the PPK method provides accuracy comparable to GCPs (Zeybek 2021).
All images were processed using the Pix4D Mapper software (Pix4D 2017).Due to cloud cover, nearly 5% of the available 1492 images could not be processed, leaving 1419 images for further processing.After the bundle adjustment, only 0.08% relative difference between initial and optimized internal camera parameters was computed.The 3D points were computed on multiple image scales (on images with half, quarter, and eighth image sizes), which is the case useful for computing 3D points in vegetation areas.More than 283 million 3D points were automatically generated, equivalent to an average density of ≈17 points/m 3 .Thereafter, the DSM and the true-orthoimage of the study area with a GSD factor of 2 (9.96 cm) were generated.We also applied automatic noise and sharpening filters, which corrected the elevation of points using the median elevation of surrounding points, and flattened only quasi-planar areas based on surface orientation.
Four representative test sites were chosen from the study area, each with varied planting characteristics (Figure 4).The test images were chosen to illustrate a variety of planting patterns of Stone Pine trees of various ages and characteristics.

Methodology
In this research, we propose a novel enhanced probability map of local maxima (P Enh L max ) that uniquely combines the newly proposed local minima (L min ) information with the previously developed local maxima (L max ) evidence in Ok and Ozdarici-Ok (2018b), which eventually supports the calculation of the orientation symmetry used for the detection of the tree stem locations.

Enhancing probability map of local maxima
In this study, L max evidence (E L max ) is collected based on a sequence of height thresholds (h i max , where In Equation ( 1), p denotes a DSM pixel at image coordinates (row, column), n is the maximum number of height thresholds utilized, and the operator "regionalmax" symbolizes the regional maxima operation (pixels connected and have a constant elevation that is greater than the elevation of their exterior boundary pixels) of the H-maxima transformation (Soille 1999) computed using the grayscale reconstruction: In Equations ( 1) and ( 2), � symbolizes the function composition, rec J; I ð Þ denotes the grayscale reconstruction which is obtained by a successive number of grayscale geodesic dilation (d) operations until no change occurs (Vincent 1993).In Equation (2), ρ where � denotes the dilation operation, D is the flat square structuring element defining 8-neighborhood connectivity, and the operator "min" returns the minimum value.In the last step, all E L max values are scaled between 0 and 1 to generate a probability map, Þ, in which brighter pixels indicate a higher probability of local maxima regions (Figure 5(b)).
To compute the L min evidence, similar to Equation (2), we also benefit from a successive number of grayscale geodesic dilations (d) performed for the grayscale reconstruction.However, at this point, we use the complement of the input DSM input for the geodesic dilations.Thereafter, we propose the new L min evidence (E L min ) of a DSM: In Equation ( 3), h j min denotes a sequence of height thresholds (j = 1, 2, •••, m) to be used in the grayscale reconstruction rec J; I ð Þ, and the operator "max" returns the maximum value after employing all predefined height thresholds (h In the last step, all E L min values are scaled between 0 and 1 to generate a probability map, Þ, in which local minima regions are indicated with brighter pixels (Figure 5(c)).
Once the probability maps for local maxima (P L max ) and local minima (P L min ) are generated, we propose to combine them to generate an enhanced version of the local maxima probability map (L Enh max ): In Equation ( 4), both P L max and P L min are already normalized between 0 and 1, therefore the output P Enh L max provides a directly enhanced probability map of local maxima regions of a DSM (Figure 5(d)).

Combining orientation symmetry with enhanced local maxima
A radial symmetry transform based on orientation is specifically developed only to identify regions with near-circular characteristics in a DSM.Therefore, (i) the enhanced probability map of a DSM's local maxima areas and (ii) orientation symmetry must be integrated to detect regions revealing both evidences.
The probability map of orientation symmetry (P OS ) combined with enhanced local maxima is computed based on a sequence of radius values (r k , where k = 1, 2, •••, K) defined in image space: In Equation ( 5), α s denotes a series of strictness parameters of orientation symmetry (α s where s = 1, 2, •••, S), and G σ is a 2D Gaussian kernel that is rotation invariant and has a fixed standard deviation (σ) applied via the convolution operation (*).The information about the center positions of circular objects in a DSM that fulfill a certain radius r k is provided in O r k , which is computed using an accumulation performed in image space: In Equation ( 6), the superscripts ({new} and {old}) denote an accumulation in image space for a radius r k .p r k indicates the unit direction vector pointing a distance specified by radius r k from p, and it is computed by Here, c is a constant vector in pixels (c = [0.5 0.5]), g ! is the unit direction vector (= ori(p) / mag(p)) calculated using the magnitude, mag(p), and orientation, ori(p) components of the gradient (Farid and Simoncelli 2004) of the input DSM.The operator ⌊.⌋ converts floating numbers to integers (i.e. the lowest preceding).A pixel p in the probability map of orientation symmetry (P OS ) reveals very high probability values if that pixel in P OS unveils both local maxima and orientation symmetry at the same time (Figure 6(b)).Therefore, we use Otsu's multilevel thresholding criterion to automatically choose the optimal threshold value (i.e.P OS > τ otsu ) to focus on regions that potentially represent the stems of trees.Based on Otsu's objective criteria, we investigate how successful the thresholds are at separating the P OS probability map into a number of classes (c q otsu where q = 1, 2, •••, C).Thereafter, we examine the success of different number of classes using the Otsu's efficiency measure one by one, and set the threshold τ otsu to the value with the maximum efficiency to generate a binary map in which isolated regions with a pixel value 1 belong to stem regions of trees, and a pixel value 0 is the background (Figure 6(c)).

Extraction of the boundaries and stems of trees
We extract the individual crown boundaries of the trees using the Chan-Vese active contour model (Chan and Vese 2001): where the input is the DSM Þ, B defines the boundary of a closed contour, and µ 1 and µ 2 are the averages of pixel values of the inner and outer boundaries of a closed contour, respectively.B l and B a represent the contour's length and the area within the contour in image space, respectively, and the inner and outer terms in Equation ( 7) are weighted by the parameters λ 1 and λ 2 , respectively.The amount of smoothness of the contour is controlled with the parameter ρ, and the contraction bias, υ, regulates the tendency of the contour to expand outwards or shrink inwards depending on the sign of the term, either positive or negative.
During the implementation, the sparse-field levelset method (Whitaker 1998) was utilized.We perform the evolution of the contour independently for each separated region, and the initialization of an active contour is automatically set from the shape of the region in which the evolution begins, limiting the evolution to a maximum number of iterations (τ iter ) if stability is not reached.Besides, image extents of each active contour to be processed are limited by generating influence regions (Ok and Ozdarici-Ok 2018a) developed to automatically merge stem regions that are closely positioned in image space (< τ dist ), thus, resulting in a more accurate final representation of trees.Finally, individual stem locations of trees are determined by searching for the greatest value of P OS within each crown boundary of trees (Figure 7).

Reference data and strategy for accuracy evaluation
The use of Hand-held Mobile Laser Scanner (HMLS) systems has become a popular technique in (medium and less closed) forest areas.A Simultaneous Localization and Mapping (SLAM)-technologyintegrated hand-held LiDAR system was used to collect the reference point cloud dataset.Before field surveying, for each test site, a total of four GNSS multi-frequency receivers (Aschtech SP 80) are deployed as georeferencing stations for HMLS point clouds.Such measurements were carried out in certain regions within each test site with the lowest forest cover to successfully solve the initial-phase uncertainty in the measurement of GNSS.In this way, signal interruption during measurements has been minimized.All measurements were supported by the continuous observation coordinate system CORS-TR technique, and the coordinates of the georeferencing points in the instant centimeter level in the field environment were determined.
Hand-held laser scanning operations were carried out on 18 August 2020, with the Geo-SLAM ZEB HORIZON system (Geo-SLAM 2021).In field surveys, while georeferencing measurements were carried out by GNSS receivers, HMLS measurements were simultaneously conducted using the GEOSLAM ZEB-Horizon SLAM system.These measurements were made at normal walking speed without any planning but by returning to the starting point by performing The object-based performance of the proposed methodology is also assessed using Equation (8) as well.If an output object has at least a 60% pixel overlap ratio with a Stone Pine tree in the reference data, it is labeled as TP (Ok and Ozdarici-Ok 2018a).If an output object does not match any of the Stone Pine trees in the reference data, it is classified as FP.If an output object matches a reference Stone Pine tree with a certain degree of overlap (i.e. less than 60%), we mark it as FN.
We also evaluate the 2D positional accuracy of stem locations detected automatically.Such positional accuracy assessment relies on a critical assumption: the greatest value of P OS within each crown boundary of trees can be compared to the actual position of the stems.Even though trees with leaning or forked stems contradict this assumption, the Stone Pine is well known for its umbrella shape, which appears to give support to that assumption.For the assessment, all reference stem locations (i.e. a center of the stem of each Stone Pine tree) are manually extracted from the HMLS reference point clouds, and subsequently used to assess our automated results.To do so, we first normalize the height of each point w.r.t.terrain height computed using the ground extractor module of the 3DReshaper software (Figure 9(b)).A distance-based segmentation step (Figure 9(c)) followed by the cylinder fitting process (Figure 9(d)) using the least-squares estimation (Figure 9(e,f)) was applied to extract precise stem locations.The Euclidean distances between the automatically extracted stem locations and the manually extracted ones are used to compute a set of measures (minimum, maximum, median, and the root mean square (RMS) distances) between the output and reference tree stems.

Results
The outcomes of the proposed tree extraction methodology are presented in this section.First, the parameters required to run our methodology are discussed.Next, numerical and visual findings are provided for the test images selected.Following that, our results are compared to the results of earlier works.Finally, the proposed framework's shortcomings are discussed.

Implementation and discussion of the user-defined parameters
MATLAB R2015 was used for the implementation of the proposed methodology.All tests were carried out using a laptop PC with a quad-core Intel i7 processor running at 2.40 GHz and 16 GB of RAM.We considerably benefit from built-in parallel processing in MATLAB during the implementation to effectively compute different parts of the methodology (e.g.all computations demanding a sequence of height thresholds, operating the active contour, etc.).
Our methodology, like any other object extraction process, needs the user to set several parameters (Table 1).The effect of each parameter on the performance measures (Precision, Recall, and F 1 -scores) was investigated in a series of experiments with scenarios for particular input parameters.
The first sequence of height thresholds, h i max where i = 1, 2, •••, n, is required for the formation of a probability map of local maxima (P L max ).For that parameter, we initiate the height sequence (h 1 max ) with the anticipated vertical accuracy of the input DSM (i.e.h 1 max ¼ 10cm), and constructed a height sequence vector with entries that are also increased by that value.The value of n was set to 8 (i.e.h 8 max ¼ 80cm) since it slightly offered the optimum balance between precision and recall measures (Figure 10(a,b)).The second sequence of height thresholds, h j min where j = 1, 2, •••, m, is required for the formation of a probability map of local minima (P L min ).In this study, by taking into account the estimated height of trees above ground level, we launch a height sequence (h 1 min ) with a value of 1 m (i.e.h 1 min ¼ 1m), and constructed a height sequence vector with entries that are also increased by that value.The value of m was set to 10 (i.e.h 10 min ¼ 10m) since it offered the best balance for the pixel-and object-based measures (Figure 10(c,d)).
The third sequence of thresholds, r k where k = 1, 2, •••, K, is required during the generation of the probability map of orientation symmetry (P OS ).We trigger the accumulation based on a radius value favored for the detection of young Stone Pine trees (i.e.r 1 = 50 cm), and break it with a radius defined according to the crown size of the frequently observed Stone Pine trees (i.e.r K = 6 m).Note also that a few numbers of trees with very large crown sizes reaching radii values around 10 m also exist in our study area.However, thanks to our accumulation process running in image space with an increment of 1 pixel (thus K = 57), we can successfully detect even such Stone Pine trees with an extremely wide radius without the requirement of increasing the value up to 10 m (Figure 10(e,f)); hence, assisting in the speeding up of processing.
The fourth sequence of thresholds, α s where s = 1, 2, •••, S, is necessary to focus successfully on radially symmetric objects.We activate the strictness value at 2 (i.e.α 1 = 2) as lower strictness values (e.g.α ≤ 3) put more emphasis on non-symmetric objects along with the radial symmetric ones.We stop the aggregation with a strictness value of 5 (i.e.α 1 = 5), with an  Default values (i.e. 1) were selected for both the weighting parameters (λ 1 and λ 2 ) of the active contour model when extracting tree boundaries.We also follow the same parameter values for the rest of the other parameters of the active contour as meticulously evaluated in Ok and Ozdarici-Ok (2018a).We adjusted the smoothness parameter ρ to 0.1, aiming for less smoothing at the tree crown borders, and therefore allowed for the collection of finer crown details.We preferred a negative signed value for the bias variable (i.e.υ) since an expansion during the evolution of the contour is required.The maximum number of iterations (τ iter ) was fixed at 200, which was found to be the threshold at which the results stabilized.Lastly, we merge the stem regions that are very closely positioned (i.e.40 cm) in object space (Figure 10(i,j)).
Figure 4 shows the total pixels in terms of row and column numbers in each test image, and Table 2 presents the elapsed time necessary for each section.The entire processing of four-test images took around 8 min to complete.The computing time of the proposed methodology is primarily influenced by two parameters.The first parameter is the total number of pixels in the input DSM, and the second parameter is the number of regions identified contributing to the final processing section.

Automated extraction
The numerical and visual findings for the test images are summarized in Table 3 and Figure 11, respectively.The results of the proposed methodology on three test images (#2-#4) appear to be robust.The pixel-based recall ratios estimated for these test sites are ≈95%, with precision ratios ranging from 82.4% to 85.3%.For test site #1, the worst pixel-based and object-based results are computed.These findings are due to two particular reasons: (i) some of the small trees do not support enough orientation symmetry and are, therefore, overlooked (marked as FN), and (ii) some of the correctly identified trees do not exist in the reference data since they are not labeled as Stone Pine trees (labeled as FP).It is worth noting that automated extraction yielded an 82% pixel-based F 1 -score ratio even with such erroneous cases.
The overall object-based precision ratio is computed to be greater than 99%, indicating that a small number of objects, mostly in test site #1, are incorrectly labeled as Stone Pine trees.This is mainly due to the Red Pine trees that also exist on the site.The overall object-based recall ratio is also successful (i.e.96.3%), with the main issue being some of the tree's small crowns, which cause inconsistencies in accumulation during the orientation symmetry.Nevertheless, the overall object-based F 1score ratio is computed to be 97.7%,considering the outcomes of all four test sites.
We must point out that the gradient orientation of a DSM that involves a coarse GSD will be less exact than the orientation component obtained from a DSM that involves a high GSD.Therefore, as long as the GSD of input DSM is comparable to the one utilized in this research, we believe our automated approach can be applied to any other geographical location consisting of Stone Pine trees.In this context, our approach requires users to set several parameters; however, an important fact is that we use a series of thresholds for a set of parameters (i.e.height thresholds (h i max , h j min ), radius (r k ), strictness (α s ) and Otsu's (c q otsu )), which provides more versatility than fixing a parameter to a particular threshold.Our sensitivity analysis carried out for different parameters proves that the results of our approach do not significantly affect from most of the parameters.However, the steepness of the terrain could be an important issue, as our approach was only tested in scenarios with relatively flat terrain.For the enhanced version of the local maxima probability map, which requires additional investigation in steep terrain, it might be possible that unexpected results could occur, and further consideration for critical parameters might be required.

2D positional accuracy
The numerical and visual findings of the automatically extracted stem locations in terms of 2D positional accuracies for all test images are summarized in Table 4 and Figure 12, respectively.We produced one-to-one mapping between the automatically extracted stem locations and the manually extracted ones after assessing the overlapping information of objects in both datasets before the calculation of the Euclidean distances.It is important to note that this process is mandatory because the    automated approaches are always prone to underor over-segmenting the trees.Therefore, the automated one-to-one mapping procedure automatically mapped 757 stems out of 964 existing in the reference data.According to the numerical findings in Table 4, the overall RMSE computed is 0.72 m in object space with a worst-case of 1.11 m for test site #1.For each test site, we see a few stem locations with relatively large Euclidean distances compared to the reference set; however, the overall median distance computed is just 36 cm in object space (less than 4 pixels), proving the success and robustness of the proposed methodology.Based on the graphical representation of errors (each vector originates from the reference stem position and extends according to the stem position extracted) provided in Figure 12, we observe almost no systematic 2D positional errors.
We observed in test site #2 that the positional errors of stems extracted for relatively small-sized trees are significantly less than those extracted for relatively tall and old trees (Figure 12).This is due to the reason that such tall trees may involve several peaks within their crown boundaries and may support multiple orientation symmetry information; thus, the probability of the precise identification of the true positions of their stems is reduced.
As well issued in recent years, multiple factors impact the quality of UAV-based outputs, point clouds, DSMs, orthomosaics, and so on (Tmušić et al. 2020).Besides, additional vegetation-related factors, e.g.vegetation leaf-off and leaf-on seasonal variation, may also influence the output accuracies (Tomaštík et al. 2019).Standardized UAV mission design parameters (e.g.flying altitude, Ground Sampling Distance (GSD), forward/side overlaps,  and in-flight calibration of sensor parameters) can be promoted without an extra effort to collect highquality datasets.Nevertheless, the need for fieldwork before UAV flight for ground control setup limits the potential benefits of timely and effective mapping.Based on our findings, we feel that a GCP-free framework using GNSS PPK for positioning appears to be successful.Nonetheless, our 2D positional accuracy is comparatively lower than previously reported PPK-related forest studies (e.g.Tomaštík et al. 2019).However, note that our result is not a direct performance assessment of a PPK-related study but an evaluation of stem locations automatically extracted after an end-to-end unsupervised object extraction framework.Therefore, our 2D positional performance eventually accumulates all  sources of uncertainties involved during the end-toend evaluation (i.e.UAV mission design, PPK quality, point cloud, and DSM generation, unsupervised framework, and reference data quality).

Comparison with other approaches
We compare our results against three state-of-the-art techniques to better understand the quality of our findings (Tables 5 and 6).It is worth noting that we used nDSMs of the test regions to meet the input criteria of those previously developed approaches, and in this study, nDSMs were successfully generated by using the approach in Pingel, Clarke, and McBride (2013).Except for the technique employed by Dalponte et al. (2015b), where we resampled the GSD of the nDSM to 0.5 m to achieve representative outputs, the GSD of the input nDSMs was precisely the same as the DSM used for our approach.We found that the default values produced the best outcomes for all techniques, and as an input set by the methods, we set the maximum radius of a tree to 10 m and the minimum treetop height to 1 m above ground if needed.
The comparative results in Tables 5 and 6 demonstrate that the proposed methodology can provide superior or comparable results in almost all of the test sites for this area.Notably, our sensitivity analysis (Figure 10) for various parameters demonstrates that the results of our approach are not significantly affected by the majority of parameters; however, we must emphasize that this issue must also be validated for different regions.When compared to the findings of the state-of-the-art techniques, the proposed method achieves overall improvements of nearly 2% and 3% for the pixel-and object-based evaluations, respectively.Among the test sites, 3 out of 4 F 1 -scores computed are found to better the results of previous studies.The recall ratios computed using Popescu and Wynne (2004) appeared to be the best for almost all test sites; however, this was achieved by sacrificing precision ratios that eventually resulted in a relatively large number of FPs, both for pixel-and object-based manner (Figure 13).
We also compare our findings to three previous methods that rely on orientation symmetry information (Tables 7 and 8).The study by Ok and Ozdarici-Ok (2018a) developed exclusively the orientation symmetry information used to extract individual trees.In Ok and Ozdarici-Ok (2017), post-processing of binary local maxima information to filter the outputs of a DSM's orientation-based radial symmetry was proposed.Later, the study by Ok and Ozdarici-Ok (2018b) combined orientation symmetry and local maxima in a probabilistic framework.Note that, for all those approaches, we used the default parameter settings except for the radii parameters (r 1 and r K ), which we set based on the radii values given in Table 1.Because of the innovative integration of orientation symmetry with the enhanced local maxima, the proposed framework fully outperforms previous approaches based on orientation symmetry, as shown in Tables 7 and 8.These results also highlight the use of new probabilistic local minima information in improving the computation of orientation symmetry.

Conclusions
A new methodology for extracting Stone Pine trees from UAV DSMs is introduced.In the first stage of the methodology, an enhanced probability map of local maxima handling the newly proposed local minima information uniquely is presented.Following that, a new probability map aiding to improve the evidence of orientation symmetry is utilized during the extraction of regions belonging to Stone Pine trees.In the final stage, the Chan-Vese active contour model, which processed each region independently, is used to recover the crown boundaries and stem locations of individual Stone Pine trees.Four test sites with distinct planting and distinct characteristics are preferred for evaluating the proposed methodology, with overall pixel-based and object-based F 1 -scores of 88.3% and 97.7%, respectively.Additionally, our findings are compared to six previously developed tree extraction approaches, and the proposed methodology's effectiveness is discussed.Additionally, the accuracy of automatically derived stem locations is compared to solid reference data collected during fieldwork in terms of 2D positional accuracies.Several stem locations with large Euclidean distances are observed during the analysis when compared to the reference data; however, the overall median distance is computed to be 36 cm in the object space (less than 4 pixels), demonstrating the proposed methodology's effectiveness and robustness.
In the future, we plan to adopt RGB images into our framework to further improve the borders of individual trees.In this sense, it would be interesting to develop a multi-source strategy that is capable of working with images having radiometric inconsistencies due to poor visibility.In this respect, we intend to expand our research to include popular learning-based techniques (e.g.deep learning) if a sufficiently large number of training datasets for Stone Pine trees can be collected.Besides, we will evaluate the effectiveness of our newly proposed methods in other regions with rugged landscapes, preferably in other parts of the Mediterranean region.The ultimate objective of this research is to create a realistic 3D representation of Stone Pine trees; therefore, we are also interested in working with innovative point-cloud processing techniques that operate with the output of HMLS systems.

Notes on contributors
Asli Ozdarici-Ok received her PhD degree from Middle East Technical University, Ankara, Turkey, in 2012.She is now an associate professor at the Academy of Land Registry, Ankara HBV University.Her research interests are in object detection, image classification for remote sensing, and applications of land and mass evaluation.
Ali Ozgun Ok received his PhD degree from Middle East Technical University, Ankara, Turkey, in 2011.He joined the Department of Geomatics Engineering at Hacettepe University in 2018, where he has been a full professor in the same department since 2021.His research interests are in the areas of photogrammetric processing and information extraction from airborne and spaceborne images, image processing, and computer vision with applications to remote sensing.
Mustafa Zeybek received his PhD in 2017.Currently, he is the head of the Department of Architecture and Urban Planning at the Güneysınır Vocational School of Selcuk University and is a faculty member.His main research areas are UAV, mobile LiDAR, terrestrial laser scanning, deformation measurements and analysis, point cloud processing, forest measurements, and road surface distress analysis.
Ayhan Atesoglu received his PhD degree at the University of Bartin, Turkey.Now, he is a professor at the Department of Forest Engineering, Bartin University.His research interests are in remote sensing and geographic information systems integrated applications, and land use/cover classification and change detection.

Figure 1 .
Figure 1.Proposed methodology for the automated extraction of individual Stone Pine trees.
investigated the changes in the Stone Pine cover for two different years in Mount Lebanon in which a hyperspectral image and a Landsat 7 image along with topographic maps were processed.Alonzo et al. (2016) tested hyperspectral imagery and waveform LiDAR data to map multiple urban forest species including Pinus pinea in downtown Santa Barbara, California.Guerra-Hernandez et al. (2016a) estimated individual Stone Pine tree parameters (location, tree height, and crown diameter) using DSM computed from a high-resolution UAV imagery in Portugal using two segmentation strategies (mixedpixel and region-based).In another study, Guerra-Hernández et al. (2016b) estimated aboveground biomass components for three types of Mediterranean trees including Stone Pine using lowdensity airborne laser scanning data in southwest Spain with two different modeling approaches.Detection of growth changes in individual Stone Pine trees was conducted by Guerra-Hernández et al. (2017) using multi-temporal UAV-derived CHM and object-based image analysis in central Portugal.Awad (2018) collected field-based spectral signatures for Stone Pine forests to test and compare different multispectral and hyperspectral satellite images for classification in two different regions in Lebanon.Demir (2018) extracted individual Pinus pinea trees by analyzing elevation flows on the CHMs derived from overlapping UAV images in a test site in the southern part of Turkey.Blázquez-Casado et al. (2019) aimed to extract two different Stone Pine tree species (Pinus pinea L. and Pinus Pineaster Ait.) at a mature development level by low-density LiDAR and Pleiades satellite images in the northern part of Spain using a random forest model.Selim, Kalaycı, and Kılçık (2020) collected height information from UAV images with the assistance of spherical astronomy in the south of Turkey for different types of trees including Pinus pinea.In a recent work, Gülci et al. (2021) evaluated the performance of UAV-based measurements for tree attributes in timber volume estimation of Stone Pine trees in the north-western part of Turkey.

Figure 2 .
Figure 2. Overview of the study area.

Figure 3 .
Figure 3. Orthomosaic of the study area, and sample images from the study area.

Figure 4 .
Figure 4. Test sites selected from the study area (The red polygon in each test site presents the coverage of the reference data collected during fieldwork; the blue triangles represent the positions of GNSS receivers deployed for HMLS measurements).

Figure 5 .
Figure 5. (a) DSM of test site #3.(b) Probability map of local maxima.(c) Probability map of local minima.(d) Enhanced probability map of local maxima.Notice how unnecessary local maxima information is suppressed in (d) compared to (b).

Figure 6 .
Figure 6.(a) DSM of test site #3; (b) Probability map of orientation symmetry combined with enhanced local maxima; (c) Output binary map.

Figure 9 .
Figure 9. (a) a sample HMLS point cloud; (b) Detection of ground points; (c) segmented objects; (d) the output of cylinder fitting process; (e) and (f) Cylinder fitted to a single stem using least squares.

Figure 10 .
Figure 10.Pixel-Based (left column) and object-based (right column) performances of the proposed methodology.

Figure 11 .
Figure 11.Visual outcomes of the proposed methodology (The results of the test sites #1-#4 are provided from top to bottom, respectively; the input DSM, the extracted regions for trees, and the visual output of the accuracy assessment are given from left to right, respectively).

Figure 12 .
Figure 12.Graphical representation of the error vectors (note the scale of error vectors shown).

Figure 13 .
Figure 13.Comparison with the previous studies (The results of the test sites #1-#4 are provided from top to bottom, respectively; the visual outputs of the approaches in Popescu and Wynne (2004), Swetnam and Falk (2014), Dalponte et al. (2015b), and ours are given from left to right, respectively; green, red, and blue colors represent TP, FP, and FN pixels, respectively).

Table 2 .
The elapsed time of each section of the proposed methodology.

Table 3 .
Pixel-and object-based performance of the proposed methodology.

Table 4 .
2D positional accuracy of automatically extracted stem locations.

Table 5 .
Pixel-Based comparison the previous approaches and the proposed methodology (the best and second-best performances across approaches are illustrated in green and orange colors, respectively).

Table 6 .
Object-Based comparison of the previous approaches and the proposed methodology (the best and second-best performances across approaches are illustrated in green and orange colors, respectively).

Table 7 .
Pixel-Based comparison of the previous approaches based on orientation symmetry and the proposed methodology (the best and second-best performances across approaches are illustrated in green and orange colors, respectively).

Table 8 .
Object-Based comparison of the previous approaches based on orientation symmetry and the proposed methodology (the best and second-best performances across approaches are illustrated in green and orange colors, respectively).