IDENTIFICATION OF PEACH TREE TRUNKS FROM LASER SCANNING DATA OBTAINED WITH SMALL UNMANNED AERIAL SYSTEM

Agricultural robotics rely on digital tools and sensor integration in order to improve efficiency and sustainability of cultivations. One part of orchard inventory is the identification of a tree trunk i.e. localization and diameter determination. However, this is a challenging task, due to thin trunks, presence of leaves and low branches. In this paper we present a case study for determining these parameters using the example of peach orchard, for which a high-density LiDAR data (over 3000 points/m) was obtained with a small unmanned aerial system (UAS) during a leafy and leafless season. We applied point thresholding by height and by components of normal vector, in order to identify points reflected from trunks. Alpha-shape algorithm was used to aggregate together points, that belong to the same trunk and their centroid determined the trunk location. Trunk diameters were calculated using two alternative approaches: the Principal Component Analysis (PCA) and circle fit. For the leafy season trunk identification is challenging. Omission errors were caused due to few reflections from trunks and commission errors occurred because of the unfiltered reflections from low branches and young twigs oriented towards the ground. All 194 trunks were identified from data collected during the leafless season. The accuracy of tree location was 0.27 m and the accuracy of diameter determination using PCA was 0.03 m.


INTRODUCTION
Modern agriculture and horticulture are becoming more and more automated, efficient and sustainable due to the use of modern technologies. Agricultural robotics, which is based on a variety of digital tool and sensor integration, is developing rapidly and results in the reduction of human resources for harvesting (Shamshiri et al., 2018). Robots and vehicles can already perform some cultivation treatment, i.e. target spraying (Ogawa et al. 2006), pruning (Ishigure et al., 2013), weeding (MacKean et al., 2017), nursing (Jørgensen et al., 2007) and general tillage (Bawden et al., 2014). Over the last decade, there has been intensive work on building an autonomous harvesting system (Li et al., 2011). Although automatic harvesting of ground crops over flat fields is already well established (Zhang et al., 2014), autonomous robotic fruit picking in orchards is still a challenging task. This is mainly due to complexity of decision and picking system, which has to identify the ripeness of fruit under different lighting, often in the presence of branches and leaves, as well as to gently collect fruit without damaging the fruit and tree (Shamshiri et al., 2018).
Identification of fruit in an orchard allows accurate estimation of the yield, knowledge of which is critical for the proper management of water resources, pesticides dosing, involvement of equipment and labour (Linker, 2017). Yield in an orchard can be roughly estimated by multiplying the average number of fruit collected from a single tree by the number of trees per unit area, i.e. orchard density (Gongal et al., 2015). Such results are reliable for large and well-structured courts (Sarron et al., 2018). Otherwise, a detailed knowledge on the exact number of trees and fruit mass for each tree is required.

* Corresponding author
In this context, some work concentrates on an accurate determination of the number and fruit mass on a single tree, using a variety of sensors. Payne et al. (2013) and Qureshi et al. (2017) counted mango fruit using daytime images of individual trees. A similar approach was successfully used for apples (Linker et al. 2012;Silwal et al., 2016;Thanh et al., 2016). Bargoti and Underwood (2017) presented an image-based object identification framework for fruit detection. Red and infrared laser diodes were used for cherries (Tanigaki et al., 2008), thermal imaging for citrus (Bulanon et al., 2008) and apple fruits (Stajnko et al., 2004). Colour (RGB) and near-infrared (NIR) imagery was captured for vineyards (Nuske et al., 2014) and sweet peppers (Sa et al., 2016). Okamoto and Lee (2009) detected green citrus from hyperspectral images. Recently, Light Detection and Ranging (LiDAR) is also used as a tool for fruit detection (Stein et al., 2016;Underwood et al, 2016).
For automated harvesting and tree determination for yield estimation it is beneficial to localize tree trunks, i.e. to determine the position of the geometric centre of a trunk and its diameter. Machine vision for detection of the trunk in orchards often fails, due to small trunk diameters, presence of branches, leaves and fruits (Gongal et al., 2015). It is particularly difficult when trees have overlapping crowns, low branches and the crown base heights are low.
Usually tree crowns are detected from remote sensing imagery (Ke, Quackenbush, 2011), but the success rate of such methods depends on tree species, crown size, distance between trunks, terrain topography and image resolution (Hirschmugl et al., 2007). Satellite multispectral imagery was used to detect citrus trees with an accuracy from 78% to 98% (Santoro et al., 2013) and oil palm trees with an accuracy of 98% (Santoso et al., 2016). Airborne imagery allowed semi-automatic detection and counting of 95% of oil palm trees (Shafri et al., 2011). Only 80% of citrus trees were detected from multi-spectral images obtained from unmanned aerial system (UAS) (Koc-San et al., 2018). Recently, a Convolutional Neural Networks (CNN) for tree identification was also applied. Satellite imagery for palm tree detection using CNN allowed to identify 92% to 98% of trees depending on a study area (Cheang et al., 2017). CNN approach combined with UAS images was also presented for citrus (Csillik et al., 2018) and orange (Osco et al., 2020) orchards.
Trees can also be detected using point clouds obtained from LiDAR. A terrestrial scanning system can map fruit distribution, detect and predict yield for an individual tree in an almond orchard . Orchard inventory using a terrestrial scanning system was also developed for an apple orchard (Gené-Mola et al., 2020). 99% of apple trees were identified and some geometric parameters were determined with a sub-decimeter accuracy using high-density point cloud obtained with UAS (Hadas et al., 2019).
In this paper a study on trunk determination in a peach orchard using a high-density point cloud obtained with LiDAR from a small UAS is presented. Results of trunk location and diameter determination are compared with direct field measurements for validation purposes.

Test area and test periods
A part of a peach (Persica vulgaris Mill.) orchard located in South West Poland near Trzebnica city (51.3100° N, 17.0457° E) was selected as the test area. The study area was about 3500 m 2 , covered by 194 trees planted in 12 rows. The trees varied in height, crown size and branch density.

Figure 1. A peach tree in the test orchard during a leafy period
Two test periods were selected: September 2017 and March 2018, to which we further refer to as leafy and leafless season, respectively. During the leafy season, the ground was covered with some dry and green grass, and weeds of different lengths ( Figure 1). The trees had no fruits but branches were full of leaves. During the leafless season only individual leaf buds on branches were present. However, under some trees parts of branches were left, that remained after pruning and the grass was much higher.

LiDAR data
We used a custom-made small UAS which consisted of hexacopter Leica Aibot X6V2, laser scanner Velodyne HDL-32E, Global Navigation Satellite System (GNSS) receiver NovAtel OEM615 and Inertial Measurement Unit (IMU) Sensonor STIM30 (Karpina et al., 2016). A GNSS receiver Leica Vica GS14 served as a ground reference station.
Data was collected during two 6 minute long flights, one for leafy and one for leafless season. For both flights, the average flight height was 30 m. UAS trajectory and georeferencing of LiDAR data was obtained from a common processing of GNSS and IMU data, following the typical procedure of kinematic LiDAR data processing described by Jozkow et al., (2017). As a result, two point clouds were obtained, one for each test period, which were georeferenced in the Polish coordinate frame PL-2000. Each point cloud consisted of more than 100 million points.

Direct field measurements
During the leafy season we visually inspected the orchard and counted the number of trees planted inside the test area. Then, we randomly selected 50 trees for which we measured trunk location. For this purpose, we used a GNSS receiver Trimble R6, which determined its position in the Real-Time Network (RTN) mode. This approach allowed us to measure the location of tree trunks with centimeter-level accuracy. Moreover, we used a ruler and measured trunk circumferences about 0.5 m above the ground. From circumferences we calculated stem diameters, which varied from 0.051 m to 0.154 m (0.114 m on average). For trees, which already had branches below 0.5 m in height, we measured the circumference just below the lowest branch.
During the leafless season we confirmed the presence of all trees in the orchard, including the 50 test trees. We assumed that they neither changed their location nor the trunk grew significantly, therefor we did not renew our measurements.

Pre-processing of point-clouds
Pre-processing of point clouds was performed in the CloudCompare v2.10.2 (www.danielgm.net/cc/) software. First, point clouds were cut to the extent of the study area, which limited the number of points to about 10 million per point cloud. Then, we used the Cloth Simulation Filter  to automatically classify points as ground and non-ground (vegetation) points. Ground point were used to generate a Digital Terrain Model (DTM) with a grid size of 0.5 x 0.5 m. Afterwards we calculated normalized heights (Hnorm), i.e. point distances to the DTM (Figure 2). Finally, normal vector components NX, NY, NZ were calculated (each in the range from -1 to +1) using plane as a local surface model, considering a spherical neighborhood with a radius of 1 meter and minimum spanning tree to define the orientation.

Strategy for trunk identification
The normalized point cloud was further processed with Matlab in order to extract points describing trunks. First, non-ground points were filtered by height. For leafless season the accepted height range was 0.2 m to 0.5 m, while for the leafy season the range was 0.1 m to 0.4 m. The higher bottom boundary of the height threshold in the leafless season was due to the presence of branches that remained after pruning, and the lower top boundary of the height threshold during the leafy season was due to many hanging branches and leaves. From two normal vector components NX, NY we took the absolute maximum for measuring the tilt of a local surface model with respect to a horizontal plane. For both test periods we removed points with the value <0.3, in order to remove reflections from nearhorizontal surfaces.
In the next step, alpha-shape algorithm (Edelsbrunner et al., 1983) with 0.2 m radius was performed on remaining points. As a result, preliminary trunk contours were obtained. The geometric centre of all points inside a contour was assumed as trunk location. Additionally, a simple quality control procedure was performed, i.e. distances between trunk locations were calculated. Due to regular planting of trees, we assumed that there should be no trunks which are closer to each other than 1 m. Therefore, if a distance between two trunks lower than 1 m was found, the trunk location was determined as the geometric center of points inside both corresponding contours.
Finally, we used points reflected from a single trunk in order to determine trunk diameter, following two approaches. In the first approach, we used Principal Component Analysis (PCA) and the diameter was the average value of ranges of the first two components. In the second approach we used a circle fit following the procedure described by Pratt (1987), and the diameter was calculated as double the radius of the fitted circle.

Validation
We calculated the total number of automatically determined trees and compared it with the true number of trees from direct orchard inspection. The accuracy of trunk location determined from LiDAR data was analysed by calculating the distances between determined trunk locations and reference field measurements for the 50 test trees. For test trees we also compared determined trunk diameters with diameters calculated from trunk circumferences, which were measured in the field. Height thresholding did not allow to filter out grass points, especially during the leafless season, when grass and weeds were relatively high. Moreover, some reflections from low branches still remained in the leafy season. We calculated inclination angle as a maximum of absolute values of the first two components of the normalized vector:

Cloud filtering and trunk location
We noticed that the inclination angle for trunk points was relatively large (Figure 3). Therefore, the removal of points with low inclination angle allowed to filter out grass points as well as some points reflected from branches in the lower part of the tree crown. Filtering results for the leafy season were not satisfactory and the corresponding trunk identification failed for many trees ( Figure  4a). We justify this by the two following factors. First of all, there were less reflections from trunks ( Figure 3a) due to leafy crowns that prevented laser beam penetration. As a consequence, omission errors were present (Figure 4a, bottom left). Secondly, low branches and young twigs oriented towards the ground were not filtered out, hence classified as trunk points, which caused obvious commission errors (Figure 4a, top right). Although we tested several other values for height and inclination threshold, there was no clear improvement of results.
For the leafless season there were sufficient reflections from tree trunks (Figure 3b) for successful identification of a trunk for all trees in the study area (Figure 4b). Identified trunk locations clearly formed rows and were separated with regular distances from each other, therefore reflecting the orchard structure well.

Validation
Locations of identified trunks were compared with corresponding location of trees measured in the field ( Figure 5). From the 50 test trees, only 29 were detected during the leafy season. The root mean square error (RMSE) of distances between identified trunks and their location measured in the field was 0.77 m, with a maximum distance of 1.26 m for tree no. 11. Only 11 trunks were identified with accuracy better than 0.50 m.
For the leafless season the RMSE of distances was equal to 0.27 m. The maximum position error was 0.66 m for tree no. 4. Only two trees were identified with accuracy worse than 0.50 m and 24 trunks of the 50 test trees were identified with accuracy better than 0.20 m. Figure 5. Distance between the measured and detected stem center location for the leafy and leafless test periods Due to the unsatisfactory results of tree trunk identification during the leafy period, validation of stem diameter determination was performed only for the leafless period. Figure  6 presents the comparison of stem diameters obtained with PCAbased approach and from circle fitting with diameters calculated from circumferences measured directly in the field for the 50 test trees. Both methods were able to determine a diameter for all test trees, but with significantly different accuracy. The RMSE of differences between diameters determined from LiDAR with circle fitting and measured in the field was 0.06 m, with the largest difference of 0.28 m for tree no. 14. For 42 trees the diameter error was below 0.05 m and for another 4 trees it was smaller than 0.10 m. The correlation coefficient between determined and measured diameters was 0.48.
Results obtained with the PCA-based method were superior to circle fitting. The RMSE was 0.03 m, maximum error reached 0.07 m and the correlation coefficient was 0.76. Locations of 45 tree trunks were determined with accuracy better than 0.05 m. Figure 6. Comparison of stem diameters from field measurements with diameters detected for the leafless period using the two methods

CONCLUSIONS
We have demonstrated that a high-density point cloud obtained from LiDAR with a custom-made UAS allows to identify the trunk of peach trees. The point-cloud processing strategy was not complex, yet effective. It was composed of the following steps: height normalization, calculation of normal vector components, point filtering, point aggregation with the alpha-shape algorithm, determination of point geometric center and PCA analysis.
For the leafless season our strategy successfully identified all 194 tree trunks in the study area. For the 50 test trees, the accuracy of tree location, by means of the RMSE of horizontal position error, was 0.27 m. Trunk diameter estimation with PCA was superior to circle fitting approach and the RMSE of trunk diameter estimation with respect to direct field measurements was 0.03 m.
A more sophisticated strategy is required for processing a noisy point-cloud obtained during the leafy season. In such a case, dense leafage of crowns prevents laser beam penetration and few points are reflected from trunks. As a result neither all trees are identified, nor is the accuracy of trunk location satisfactory. Classification of trunk points in the presence of various disturbances in a complex and heterogeneous environment (e.g. high grass and nettle, low branches, young twigs oriented towards the ground) still remains a challenge.