Canopy Parameter Estimation of Citrus grandis var. Longanyou Based on LiDAR 3D Point Clouds

Abstract

The characteristic parameters of Citrus grandis var. Longanyou canopies are important when measuring yield and spraying pesticides. However, the feasibility of the canopy reconstruction method based on point clouds has not been confirmed with these canopies. Therefore, LiDAR point cloud data for C. grandis var. Longanyou were obtained to facilitate the management of groves of this species. Then, a cloth simulation filter and European clustering algorithm were used to realize individual canopy extraction. After calculating canopy height and width, canopy reconstruction and volume calculation were realized using six approaches: by a manual method and using five algorithms based on point clouds (convex hull, CH; convex hull by slices; voxel-based, VB; alpha-shape, AS; alpha-shape by slices, ASBS). ASBS is an innovative algorithm that combines AS with slices optimization, and can best approximate the actual canopy shape. Moreover, the CH algorithm had the shortest run time, and the R² values of V_CH, V_VB, V_AS, and V_ASBS algorithms were above 0.87. The volume with the highest accuracy was obtained from the ASBS algorithm, and the CH algorithm had the shortest computation time. In addition, a theoretical but preliminarily system suitable for the calculation of the canopy volume of C. grandis var. Longanyou was developed, which provides a theoretical reference for the efficient and accurate realization of future functional modules such as accurate plant protection, orchard obstacle avoidance, and biomass estimation.

1. Introduction

Citrus grandis var. Longanyou is distributed in Guang’an city, Sichuan Province, from the hills of central Sichuan to Huaying Mountain in the east. It is one of the national agricultural products of geographical indication in China; tree shapes include open-hearted, round-headed, and double-layered and Y-shaped [1]. The canopy adsorbs light, interacts with the external environment, and is the main place where fruit trees perform respiration and photosynthesis. Information perception studies have shown that the crown determines the growth and development of trees, which in turn affects their yield and economic benefits [2]. Canopy characteristics have important guiding significance for research related to trees, because they can usually be used for growth assessment [3], biomass estimation [4,5], pruning effect analysis [6], yield prediction [7], health assessment early disease detection [8,9] and in the use of plant protection products such as insecticides [10]. Fruit tree canopy characteristics include canopy height, width, and volume. Among these, canopy volume has an important influence on the nutrient demand and fruit yield of fruit trees [11]. In addition, accurate canopy volume information provides basic data for decisions related to variable spray application, and also affects the precise and variable application of pesticides [12]. Therefore, calculating the geometric parameters of the canopy quickly and accurately is of great significance for the intelligent and accurate management of orchards and plantations.

When precise information about the characteristics of the canopy is acquired, these data can be used for the fixed-point management of crops, thus improving crop yield [13]. At present, methods used to measure canopy volume mainly include manual measurement and various non-contact automatic measurement methods. Manual measurement usually involves collecting canopy height and diameter data for individual fruit trees on site; then, the canopy is treated as a regular solid shape based on the shape of the canopy, and canopy volume is calculated using a formula of approximate regular geometric volume [13]. Because canopy shapes vary from tree to tree, it is difficult to apply a standard uniform and regular geometrical shape to calculate canopy volume, and many researchers have independently proposed standards. Schinor et al. [14] assumed that a fruit tree takes on an ellipsoid shape to calculate the canopy volume. The ellipsoid model is also a common and frequently used geometric model [15,16]. In addition, managers of most Brazilian groves often calculate canopy volume by calculating the volume of a cube surrounding the canopy as the required volume; Scapin adopted this model when studying citrus canopies, independently [17]. Another type has been proposed by Li et al. [18] to measure the circumference at different heights of fruit tree canopies using a tape ruler, and then calculating the area of the circle. The volume of the fruit tree canopy is calculated by the layered volume of the cylinder measured at various heights, and the total volume is the sum of the volume of each layer. This method is more accurate and can reflect the volume of actual fruit trees. Manually measurement using the contact method is simple, easily undertaken, and requires limited mathematical knowledge for growers, but is time-consuming, labor-intensive, and inefficient; in addition, it is easy to damage a tree. Moreover, the accuracy of measurement is susceptible to mismeasurement of subjective factors. The lack of branches in the structure parameters of the system description fails to meet the requirements for an orchard measurement operation [13,19].

Approaches based on photoelectric technology have become the mainstream method used to establish three-dimensional (3D) models of objects without direct contact; this method is now used to study the canopy parameters of fruit trees. Ni et al. [20] reconstructed the canopy structure of Croton, pepper, and lemon using a monocular stereo vision system, whereas Hobart et al. [21] realized the 3D reconstruction of an apple orchard map with a RGB camera and calculated the tree trunk height. Dong et al. [22] used a high-resolution camera to achieve 3D reconstruction of apple and pear trees, and to calculate tree canopy area and diameter. Jurado et al. [23] reconstructed olive tree morphology and estimated plant height and volume using voxel-based methods by incorporating multi-spectral sensors based on high-resolution cameras. The advantages of using cameras and other devices to obtain data are their low cost, easy system development, and access to texture information. The disadvantages are that a camera is easily affected by natural light and post-processing is complicated. The advantages of a depth camera are wide-angle viewing and strong anti-interference ability, which makes up for the defects of cameras and other equipment. Li et al. [24] and Vazquez-Arellano et al. [25] used a depth camera to establish a 3D model of maize plants, and the latter applied it to plant positioning. Zhou et al. [26] realized the 3D reconstruction of the trunks, branches, and canopy of trees indoors and outdoors using a Kinect V2 depth camera, with a method featuring rapid timing and high accuracy. However, the main disadvantage of depth cameras is the low resolution of depth images. As a precision measuring approach, light detection and ranging (LiDAR) has the advantages of rapid image acquisition of 3D point cloud information and high precision, and is not easily affected by the external environment. As a result, it is currently widely used in agricultural production. Medeiros et al. [27] obtained point cloud data with a 2D laser sensor to reconstruct a model of apple trees, and then constructed the skeleton model of fruit trees for pruning. Qiu et al. [28] constructed a 3D map of cornfields and calculated the row spacing and plant height using LiDAR combined with the Global Navigation Satellite System on a mobile platform. Yun et al. [29] selected portable mobile LiDAR to obtain oak data to document hurricane and cold damage in three different sites, and evaluated the sensitivity of different rubber trees to damage from natural disturbance by calculating the included angle between the trunk and primary branches, tree height and canopy volume, and other biological characteristics.

Two methods have mainly been used to calculate canopy volume based on point clouds obtained by LiDAR: The first method involves using triangulation of points outside a point cloud to represent the object surface. This method allows the researcher to form an outer envelope polygon from the outside to the inside of the crown surface, and then calculate the volume. The second method is based on discretization to create small regular geometric or patch models within a point cloud structure [30]. The first method can be divided into two related methods, called convex hull and alpha shape. Xiao et al. [31] used a convex hull algorithm to calculate the volume of sugar beet for predicting biomass, and Lin et al. [32] proposed an improved convex hull algorithm to reconstruct the canopy of 20 trees and calculate the projected area and volume of the canopy, which greatly reduced the time required by the algorithm. Colaço et al. [30] used an alpha-shape algorithm to calculate the canopy height and volume of a commercial citrus orchard. Vauhkonen et al. discussed the application of the alpha-shape algorithm to estimating forest stock volume [33] and proposed an improved alpha-shape algorithm to estimate forest biomass [34]. Wang et al. [35] undertook canopy reconstruction and volume calculation for trees in orchards using the same algorithm. Yao et al. [36] used 2D convex hull and alpha-shape algorithms to calculate forest canopy area and volume, respectively, for the identification of forest tree species. In addition, 3D convex hull and alpha-shape algorithms were used to calculate forest canopy volume for comparison with field-measured values and the ellipsoid canopy volume in the study of Korhonen et al. [15]. The result showed that the estimated value based on LiDAR is usually smaller than the measured value. Jiang et al. [3] also calculated the canopy volume of blueberry bushes for yield assessment using these two algorithms.

Slicing- and voxel-based algorithms are based on discretization. The former borrows concepts from calculus along the horizontal Z-axis, which is split into slices; in turn, the subdivision volume is obtained using other methods, and finally the accumulated canopy volume can be obtained. The latter involves dividing the point cloud of trees into small cubes to reconstruct the trees’ canopy and then calculating the canopy volume. Wu et al. [37] used a voxel-based algorithm to reconstruct street trees, and calculated tree height, crown width, and other parameters. Zhang et al. [38] established a 13A high-spindle-shaped Malus domestica Borkh. digital model. Moreover, a cubic unit segmentation of the model was performed and the light capture rate of various branches in the total canopy and the canopy itself were analyzed. Hosoi et al. [39] established an accurate tree body model based on voxels for Camellia sasanqua and Deutzia crenata, and studied estimation methods used to accurately calculate leaf area density (LAD) and the cumulative leaf area index (LAI) under different conditions. Many researchers integrated the use of a variety of methods to reconstruct the canopy. Olsoy et al. [40] calculated the volume of sagebrush (Artemisia tridentata) by convex hull and voxel algorithms to estimate total and green biomass. Fernández-Sarría et al. [41] used convex hull, convex hull by slices, and voxel-based algorithms for street trees to reconstruct a canopy model and calculate canopy volume. Bi et al. [42] used the same algorithms for the Amorpha canopy to calculate canopy volume and verified that a voxel-based algorithm was most suitable. Cheein et al. [43] calculated the volume of the olive canopy using the same three algorithms compared with a cylindrical fitting model, and analyzed the characteristics of the different methods. Yan et al. [44] proposed another self-adaptive slicing concave-hull algorithm and compared the results with those of a convex hull, convex hull by slices, alpha-shape, voxel-based, and geometric models, which realized the reconstruction and volume calculation of a street tree canopy.

As a national agricultural product of geographical indication in China, published research studies about estimation of the canopy characteristics, and different calculated volume values using different algorithms, of C. grandis var. Longanyou have rarely been reported. In addition, few studies have been published to date on the method of calculating the crown volume by combining the alpha-shape algorithm and slices (alpha-shape by slices), or about comparative analysis based on the quantitative relationship of the calculated volume combined with the running time. In this study, based on 3D point cloud data collected by LiDAR, canopy point clouds of fruit trees were obtained by filtering and single tree segmentation. The existing manual method (MM) and five canopy reconstruction algorithms based on point clouds (APCs, specifically convex hull, CH; convex hull by slices, CHBS; voxel-based, VB; alpha-shape, AS; alpha-shape by slices, ASBS) were compared and analyzed by calculating results and determining the run time required for analysis. The relationship between different volume values is explored, and suitable algorithms are proposed for different job requirements. We calculated crown volume by combining the alpha-shape algorithm and slices (alpha-shape by slices, ASBS), and in this paper discuss the relationship among different algorithms. It was found that ASBS can better remove internal holes and gaps, and more closely approach the real canopy shape. Furthermore, V_ASBS can be calculated indirectly using a CH algorithm with the shortest running time via good quantity relationships. In addition, suitable algorithms are proposed for different job requirements. The six calculation methods shown in this work provide canopy-related geometry in addition to theoretical references designed to support intelligent orchard management in the future. The research method adopted in this study is described in Section 2. The discussion of Section 3 combines the calculation results and algorithm principles to conduct an in-depth analysis, discussing the reasons for the results and the main findings of this research. Section 4 summarizes the relationship among the volume values of different algorithms of C. grandis var. Longanyou and the optimal algorithms for different work requirements; it also emphasizes the potential contribution of this study and the improvements needed in future research.

2. Materials and Methods

2.1. Test Area and Test Equipment

The test site was a grapefruit plantation in Guang’an city, Sichuan Province, China. Citrus grandis var. Longanyou trees were used with row spacing of 5–6 m, plant spacing of about 6 m, and a tree height of about 4 m. The branches at the top are of different heights and diverge to both sides, and the branches of adjacent trees intersect each other inside, resulting in numerous holes and gaps between branches and leaves. The open-hearted shaped trees were cultivated using the ridge-tillage cultivation method; the ridge spacing is 5–6 m and the ridge height is 0.4–0.5 m. The experiment was conducted in Guang’an township, Guang’an district, Guang’an city in 2019. The local (30°50′ N, 106°76′ E) average elevation is 320 m, and the annual average temperature is 17.58 °C. The test site is in a typical hilly area of southwest China with high temperatures, a rainy climate, and medium soil fertility. A total of 36 grapefruit trees in five rows were selected as test samples in this experiment; the test scene is shown in Figure 1.

The hand-held 3D LiDAR scanning system consisted of a laser (VLP-16 LiDAR, Velodyne Lidar Co., San Jose, CA, USA), an inertial measurement unit (MTI-30-2A8G4, Xsens Co., Enschede, Holland), and a microcomputer (BRIX GB-BRi7H-8550, GIGABYTE Co., Zhejiang, China). Both are powered by a 12 V Guang Yi Sheng power lithium battery with a capacity of 30 Ah, and the LiDAR system can work continuously for about 6 hours with the power supply. According to the study by Jaakkola et al., the ranging error of a VLP-16 LiDAR scanner is 3 cm, which can meet the needs of forestry applications, and this study used manufacturing values without calibration [45]. The device connection diagram and specific parameters are described in detail in Zhou et al. [46]. The experimenters used a handheld laser radar scanning device to scan the orchard by walking a closed loop path and constructing a 3D map of the orchard (Figure 2).

2.2. Orchard Point Cloud Data Processing

The flow of processing the orchard point cloud data is shown in Figure 3. First, a 3D map of the orchard was constructed using the LiDAR scanning system; next, open-source software, CloudCompare™ (vers. 2.9.1, General Public License software, http://www.cloudcompare.org/, 9 May 2021, CC), was used to segment the test area. To reduce the amount of data processing, the ground point cloud was filtered by the cloth simulation filter (CSF) algorithm [47]. A Euclidean clustering algorithm was then used to segment a single fruit tree, which was convenient for subsequent processing. Each method, MM and APCs (CH, CHBS, VB, AS, and ASBS), was used to realize C. grandis var. Longanyou canopy 3D reconstruction and calculate each tree’s volume. The relationship among the volume values of different algorithms and their adaptation to different scenarios were explored by analyzing different methods of calculation results and the run time.

2.2.1. Ground Filtering

The 3D original point cloud map of the orchard contained both test and non-test areas. CloudCompare™ software was used to select the test area (Figure 4a). To reduce the calculation time for the subsequent processing of point cloud data, the ground point cloud was filtered by appropriate methods based on the ground terrain features. The most commonly used algorithms include the random sample consensus (RANSAC) and CSF algorithms. The RANSAC algorithm is widely used in point cloud plane fitting and segmentation, which is suitable for flat and gentle terrain; the CSF algorithm is suitable for uneven and steep terrain. Because the C. grandis var. Longanyou plantation is located in a hilly area, the ground has an uneven slope; grapefruit is typically grown using the ridge tillage method; the red wavy dotted line in the left part in Figure 4a is a topographic map of the actual ground. Because the terrain undulates, this study used the CSF algorithm to filter the ground point cloud.

The CSF algorithm in CC software has three parameters: Cloth resolution refers to the mesh size of the virtual cloth used to cover the terrain. The smaller the value, the finer the digital terrain model. Max Iterations is the maximum number of iterations of the terrain simulation. The higher this value, the higher the accuracy, but a higher Max Iterations value will require a longer time for calculation. Classification threshold refers to the threshold that divides the point cloud into ground and non-ground parts based on the distance between each point and the simulated terrain. The smaller this value, the smaller the number of ground points filtered out. The three parameter values set in this study were 0.5, 500, and 0.5, respectively. The effect after filtering is shown in Figure 4b. Most of the ground point clouds were filtered out, although there were a few sparsely scattered ground point clouds (inside the red ellipses around dotted lines); this did not affect the extraction of individual fruit trees.

2.2.2. Individual Tree Segmentation

After filtering the ground point clouds in the test area, an appropriate clustering algorithm should be used for point cloud segmentation in combination with crop growth characteristics to obtain a point cloud data of each individual fruit tree. Clustering algorithms based on 3D point clouds mainly include the Density-Based Spatial Clustering of Application with Noise [48], K-means clustering [49], and European clustering [50]. European clustering was used to extract individual fruit trees because the C. grandis var. Longanyou trees had been planted with uniform spacing and had little shade between each other. The method of searching for the nearest neighbor based on the Euclidean distance uses the kd-tree data structure. The decisive parameters in a Euclidean clustering algorithm are MinClusterSize and MaxClusterSize, which are the minimum and maximum number of points in the cluster (or in this case, representing a fruit tree), respectively. The smaller these two values, the more types are included after segmentation. ClusterTolerance refers to the tolerance variable of the distance between the points in the cluster. If the tolerance value is too small, the label marked as one tree can be considered to represent multiple trees. Conversely, if the tolerance value is too large, multiple trees will be marked as a single tree. According to the actual planting tree-row spacing and plant spacing of C. grandis var. Longanyou within each row, the three parameter values in this article were set as 50,000, 5,000,000, and 0.05 m, respectively. The clustering results shown in Figure 5 show each tree represented by a different color. Simultaneously, the European clustering algorithm also removes any residual ground points from the cloud (Figure 4b).

2.3. Calculation of Canopy Height and Crown Width of Fruit Trees

The crown height (H) and crown width (W) of fruit trees are both important characteristic variables of the crown structure (Figure 6), which directly affects the productivity and vitality of fruit trees, and can also reflect the yield of fruit. Obtaining crown height and width is convenient for orchard management, and in some countries managers determine pesticide dose according to crown height [13]. Forest canopy height refers to the distance from the highest point of the canopy to the first live branch above the ground. For C. grandis var. Longanyou, because its tree height is much lower than tree height in natural forests and the branches droop to touch the ground, the crown height of grapefruit was recorded as the distance from the highest point to the lowest point of the fruit tree point cloud according to Yang’s definition [2]. Therefore, crown height was calculated using Equation (1):

$$H={\mathrm{Z}}_{\max }-Z_{\min }$$

(1)

where Z_max and Z_min are the maximum and minimum values of Z-axis point cloud, respectively.

Crown width refers to the average north–south and east–west widths of trees, and is usually used to indicate the specifications of trees and fruit tree [51]. This study used Equation (2) to calculate canopy width:

$$W=\frac{D_{NS}+D_{EW}}{2}$$

(2)

where D_NS is the north–south width and D_EW is the east–west width.

2.4. Canopy Reconstruction and Volume Calculation

In this study, MM and APCs (CH, CHBS, VB, AS, and ASBS) were used to reconstruct the canopy of fruit trees, with the aim of calculating the canopy volume. All programs of APCs were run in MATLAB 2016A (Mathworks, Natick, MA, USA).

2.4.1. MM Method

Because C. grandis var. Longanyou are a type of citrus tree and because many researchers have often chosen an ellipsoid as the geometric model of a typical citrus tree canopy [14], this study used an ellipsoid model to represent the canopies of C. grandis var. Longanyou trees. The model diagram (Figure 6) used a blue ellipsoid for the fitting geometry of this manual measurement model, and the volume was calculated using Equation (3):

$$V=\frac{\pi \cdot W^2\cdot H}{6}$$

(3)

where H and W were measured three times using CC software to calculate the average volume [52].

Figure 7a provides a schematic diagram of canopy measurements.

2.4.2. CH Algorithm

Researchers have used four main methods to calculate tree canopy volume using the convex hull method: incremental, gift wrap, divide and conquer, and the quickhull method. The quickhull method was used in the present study [53]. This algorithm removes any points that are not within the bounds of each closed convex hull (that is, any internal points), and represents each tree crown as a 3D convex polyhedron with a triangular surface (Figure 7b). The volume determined by a CH algorithm contains all of the simplexes belonging to the smallest convex set containing point data.

2.4.3. CHBS Algorithm

The method using a CHBS algorithm mainly segments a point cloud in the Z direction based on a certain number (N) of slices; then, the point cloud is divided into N parts, and the CH algorithm is performed for each segmented point cloud. Finally, the sum of the convex hull volumes of all slices cumulatively is considered the volume of each canopy. Different slice numbers influence the experimental results differently. Because at least three points are required for each slice when calculating the volume with a CH or AS algorithm, Yan initially chose 15 cm spacing in the study of slices [44]. Furthermore, $\overline{H}$ is 3.302 m in this study (

Table 1. Comparative analysis of crown height obtained using MM and APCs.

H	Manual (m)	Algorithm (m)	Absolute Error (m)	Relative Error
Mean	3.302	3.432	0.167	0.054
Var.	0.146	0.099	0.019	0.002
Max	4.080	4.262	0.575	0.255
Min	2.250	2.825	0.002	0.001

). Therefore, to ensure that the two algorithms cover more points on the slice and facilitate the following comparative analysis, the N value in this study was 20. Figure 7c shows the effect diagram of crown reconstruction using a CHBS algorithm. The influence of different slice numbers on the calculation of canopy volume and operation performance of the CH algorithm is discussed in detail in Section 3.2.5 below.

2.4.4. VB Algorithm

The VB algorithm produces a representation of the point cloud with a small cube having certain specifications; the algorithm counts the number of cubes containing the canopy point cloud. Given the volume of the small cube, the volume of the canopy can be calculated by the number of cubes and the volume of each small cube [39]. For any irregular tree crown, a VB algorithm only needs to consider the size and validity of each cube, without considering the crown shape. Wei [54] considered that the canopy had many gaps rather than being a single entity, and used the VB method to calculate the canopy volume. The results show that the calculated canopy volume tends to be stable when the voxel size is equal to one-tenth of the canopy width. Because $\overline{W}$ in this study is 3.302 m (

Table 2. Comparative analysis of crown diameters obtained using MM and APCs.

W	Manual (m)	Algorithm (m)	Absolute Error (m)	Relative Error
Mean	3.336	3.348	0.105	0.031
Var.	0.088	0.114	0.021	0.001
Max	3.825	4.340	0.769	0.215
Min	2.729	2.837	0.004	0.001

), Lecigne [55] confirmed that a smaller voxel size is more suitable for observing subtle changes in spatial exploration than a larger voxel size in the results obtained with simulation data. Based on the specific parameters of C. grandis var. Longanyou, the voxel size of this study was 0.1 m; the VB algorithm effect diagram is shown in Figure 7d.

2.4.5. AS Algorithm

The AS algorithm is a concave hull algorithm. In this method, the simplex of the underlying triangulation is compared with the specified α, deleting the simplex with an empty external sphere and a square radius greater than the defined α. Then the volume of the 3D object is calculated [56]. It should be noted that when α is large enough, the canopy reconstruction effect will be similar to that of a CH algorithm.

The parameter α is used to determine the level of detail of the obtained triangulation, which in turn controls the fineness of the generated polyhedron, which directly affects the surface reconstruction results [15,57]. Because the slender branches of tree number 1 are more scattered than those of other trees, its crown contour contains a large amount of detailed information, which is more representative than other trees. Figure 8 takes tree number 1 as an example, showing the canopy reconstruction effects under different values for α. When α = 0.05, the canopy structure of fruit trees is expressed in a detailed manner, which expresses the canopy topological structure and degree of density. When α is 0.25 and 0.5, smaller pores in the fruit tree canopy are surrounded by triangles, which only show the spatial extension of the canopy structure, but not the detailed branching characteristics. When α = 0.75, the canopy is completely wrapped, generating a convex hull of the external point cloud, and is unable to express the information in the internal structure. According to the above analysis, by setting different values of α, different models of a fruit tree canopy can be obtained, and different 3D structures can be reconstructed, which can reflect different spatial morphologies, degrees of density, and branching characteristics of the canopy. By comparing the effect diagrams of different values and CH algorithms, it can be clearly found that when α = 0.75 (i.e., far larger than 0.05) the canopy reconstruction diagram is similar to the diagram of the CH algorithm. This phenomenon is also the same as the principle of the two algorithms. In this study, α was selected as 0.05; the effect diagram of the AS algorithm is shown in Figure 7e.

2.4.6. ASBS Algorithm

The principle of the ASBS algorithm is similar to that of the CHBS algorithm. The only difference is that the AS algorithm was used for canopy reconstruction and volume calculation for the sub-parts equally divided along the Z-axis. As with the CHBS algorithm, the value of N was set at 20 in this study; the effect diagram of the ASBS algorithm is shown in Figure 7f. The number of slices in the ASBS algorithm is also discussed in detail in Section 3.2.5 below.

3. Results and Discussion

3.1. Crown Height and Crown Width

The mean, maximum, and minimum, in addition to the variance, absolute error (AE), and relative error (RE) of MM and APCs for the crown height (H) and width (W) of 36 C. grandis var. Longanyou trees are shown in Table 1 and Table 2.

Table 1 and Table 2 show that the RE difference between the MM and APCs results was very small, and the RE of W was even smaller (the RE of H was 0.054, and the RE of W was 0.031). This indicates that the collected point cloud data are similar to the actual data, and the method used can replace the MM to a certain extent, providing a more accurate data basis and a reliable method for the calculation of canopy volume below.

3.2. Crown Volume

The canopy volumes of 36 C. grandis var. Longanyou trees were calculated using MM and five different APCs (Figure 9). The purpose of this step was not to verify the volume estimation, but to compare different canopy reconstruction methods, combined with the algorithm principle, to understand the differences resulting from different algorithms applied to the canopy of C. grandis var. Longanyou. Then the quantitative relationships of volume values among different algorithms, and the selection of appropriate algorithms for different operation requirements, is discussed.

3.2.1. Comparative Analysis of Six Different Methods

The relationship of measured volumes for a single plant was V_MM > V_CH > V_CHBS > V_VB > V_AS > V_ASBS (Figure 9). Among these, the calculated canopy volumes from the MM model were greater than those of all canopy volume calculations by APCs (${\overline{V}}_{\mathrm{MM}}$ = 19.60 m³); V_VB was about half of V_CH (${\overline{V}}_{\mathrm{CH}}$ = 17.45 m³, ${\overline{V}}_{\mathrm{VB}}$ = 8.25 m³); V_ASBS was about one-third of V_VB (${\overline{V}}_{ASBS}$ = 2.74 m³). The significant differences and algorithm characteristics of these calculation results are closely related to the structural characteristics of the canopy of a C. grandis var. Longanyou tree.

MM treats the tree crown as an ellipsoid and then the canopy volume can be calculated. In addition to the internal holes in an actual canopy, this method also includes the external space between the canopy and the theoretical model (Figure 6). Of particular note, it appears that there are fewer internal holes and gaps in Figure 6 because the canopy of C. grandis var. Longanyou is only seen from one angle. In fact, the internal holes and gaps in the canopy are much larger than those shown in the figure. Therefore, the volume value calculated by MM method is much larger than the real volume of the canopy. The CH algorithm establishes a minimum convex hull based on the canopy point cloud as a whole (Figure 7b), and holes and gaps between branches and leaves are regarded as part of the canopy; this method also creates a larger calculated volume than the actual volume of the canopy. The CHBS algorithm divides the entire canopy into several sub-parts, and establishes the sub-part of the convex hull with partial point clouds as the calculation object; however, the algorithm still contains many holes and gaps, and the calculated volume is still larger than the actual volume of the canopy. Because the VB algorithm treats the point cloud of branches and leaves as a carrier to build small cubes, the internal holes hidden by external branches and leaves can be removed when using this algorithm; as a result, the calculated volume is greatly reduced. However, because 3D cubes are formed, the final volume will be larger than the result of an AS-related algorithm. The AS algorithm can form a smaller envelope network with a reasonable setting of α, and most holes and gaps are consequently removed to form the actual contours of the canopy. Therefore, V_AS is smaller and V_ASBS is minimized.

The mean and variance of canopy volume calculated by MM are the largest, whereas those calculated by the ASBS algorithm are the smallest. In addition, the variance of the six methods decreases and tends to be stable as the volume decreases, which indicates that the holes and gaps between branches and leaves are the main factors affecting the canopy volume.

The volumes calculated by the six methods were fitted for analysis to explore the relationships between the different methods (

Table 3. Equations relating crown volumes calculated by six different methods and R².

Equation/R2	VMM	VCH	VCHBS	VVB	VAS	VASBS
VMM		y = 0.954x + 2.965	y = 1.076x + 5.083	y = 2.069x + 2.544	y = 4.891x + 3.027	y = 5.864x + 3.536
VCH	0.699		y = 1.186x + 1.449	y = 2.382x − 2.196	y = 5.492x − 1.163	y = 6.392x − 0.066
VCHBS	0.614	0.975		y = 2.008x − 3.066	y = 4.586x − 2.045	y = 5.295x − 1.013
VVB	0.515	0.898	0.920		y = 2.410x + 0.090	y = 2.773x + 0.649
VAS	0.432	0.719	0.722	0.876		y = 1.151x + 0.236
VASBS	0.470	0.733	0.724	0.874	0.996

Note: manual method, MM. Canopy reconstruction algorithms include: convex hull, CH; convex hull by slices, CHBS; voxel-based, VB; alpha-shape, AS; alpha-shape by slices, ASBS.

). The R² of MM and APCs were relatively small, between 0.432 and 0.699. The R² between the five canopy-reconstruction and volume-calculation algorithms based on point clouds were between 0.719 and 0.996. The regression coefficients of each algorithm in Table 3 are consistent with the numerical relationships in Figure 9. Next, a comparative analysis among the APCs used here is described in detail.

3.2.2. Comparison between the MM and APCs

The regression coefficients of V_MM and V_CH are close to 1, indicating that the two values are very close, but their R² is small (only 0.699). The is because an ellipsoid entity was used to simulate the canopy volume, which contains not only the holes and gaps inside the canopy, but also includes a small amount of space outside each canopy, resulting in excessive results and a large error in the estimation of canopy volume [18]. The increasing regression coefficient of V_MM and other volume measurements means that the difference between V_MM and other values increases, and the R² of linear fitting between V_MM and other volume measurements decreases with an increase in the difference in volume. This occurs because the ellipsoid model and CH algorithm overestimates the canopy volume compared with the actual value, and the canopy contains a large number of holes, whereas the VB and AS-related algorithms reduce a larger number of holes, and therefore the difference between the calculated volume of MM and APCs steadily increases, which leads to the above-mentioned trend of the regression coefficient and R². When Li et al. [18] studied the volume of the canopy of citrus trees, they found that the envelope volume measured by the most external branches and leaves of the canopy was measured by MM, whereas the canopy volume calculated accurately by APCs through the spatial distribution of branches and leaves was closer to the actual volume. Therefore, the canopy volume calculated by APCs was generally smaller than V_MM, which is consistent with our research results.

In summary, the MM model requires measures of the parameters for each tree individual and the selection of an appropriate geometry model. The differences in tree canopies make it difficult to have a unified geometric model suitable for replacing the actual canopy shape, and external space will exist between any geometric model and the real shape of the canopy, which will be much larger than the real volume of the canopy. Therefore, the calculated MM volume is far less accurate than those found using the five APCs.

3.2.3. Comparison of Canopy Volume Measured by APCs

Although V_VB is approximately 47.3% and 61.1% of V_CH and V_CHBS, respectively, the R² values are 0.898 and 0.920, respectively. The essence of the CH-related algorithm is to create an entity surrounded by triangles whose vertex is the outermost point of the point cloud data obtained by LiDAR. This volume includes all of the internal holes and the outer space between the branches and leaves; as a result, it overestimates the apparent volume. Even if the CHBS algorithm is used, the holes between the slices can be eliminated, but not effectively removed. In the VB algorithm, when no point cloud exists in these cubes, the voxels are not considered when calculating the total volume, which is why V_CHBS is larger than V_VB.

In the present study, V_AS and V_ASBS were approximately 41.1% and 33.2% of V_VB, respectively, and the R² values were 0.876 and 0.874 when comparing V_VB to V_AS or V_ASBS, respectively. Because the AS algorithm forms 3D models of all point clouds surrounding the canopy from the outside to the inside, each vertex is composed of point clouds, whereas the VB algorithm has a fixed-size voxel formed with the point cloud as the center. Therefore, the spatial volume of the voxelized 3D canopy model is larger than the polyhedron formed by the AS algorithm. As V_ASBS becomes smaller, the above volume relationship exists and R² decreases slightly (from 0.875 to 0.874).

The volume of the CH-related algorithm was much higher than that of the AS-related algorithm, and the regression coefficient between them was approximately 5; R² was around 0.72. Although both the CH and AS algorithms use a 3D model to form the canopy from the outside to the inside, the CH algorithm uses a simplex of minimum convex triangulation, whereas the AS algorithm is a concave shell algorithm that uses the parameter α to limit the number of simple bodies in the triangulation to obtain more detailed shapes [15]. Therefore, when forming the envelope network on the surface of the canopy, the AS algorithm can effectively remove the holes and gaps inside the canopy, forming a more appropriate 3D model that is more consistent with the actual shape of the canopy. Based on the above discussion, the regression coefficients of the CH-related and AS algorithms are far larger than 1. In addition, because V_ASBS is the smallest with the addition of slices, the regression coefficient of the CH and ASBS algorithms (6.392) is greater than that between the CH and AS algorithm (5.492).

3.2.4. Running Time Analysis of APCs

Five APC programs were run in MATLAB 2016a software on a Lenovo computer (processor Intel Core™ i7-8550U @1.80 GHz, 16 GB DDR4 2666 MHz memory, 256 GB SSD), and all other applications were shut down while running the program. To better compare and analyze these five algorithms, the run times of different algorithms were also calculated (Figure 10).

Figure 10 clearly shows the relationships between the run times for the five algorithms. These run times were ${\overline{T}}_{\mathrm{CH}}$ = 0.277, ${\overline{T}}_{\mathrm{AS}}$ = 1.699, ${\overline{T}}_{CHBS}$ = 2.343, ${\overline{T}}_{ASBS}$ = 3.909, ${\overline{T}}_{\mathrm{VB}}$ = 4.867. Because CH and AS algorithms essentially form a 2D plane successively, whereas the VB algorithm forms a 3D voxel, the run time of the VB algorithm was the longest. The order of stability from strong to weak was CH, CHBS, AS, ASBS, VB. Because the number of point clouds in each tree was different, the number of small cubes generated by voxel was significantly different for point clouds with different structures, whereas the differences in 2D planes generated by the CH and AS algorithms were very small. Therefore, the stability of the VB algorithm was minimal; specifically, the variance was the largest (0.974). The slice variance increases because the point cloud distribution of different layers varied and the run time for forming the plane was disparate.

Based on the reconstruction effect, volume, and run time, the MM volumes were far higher than the actual volumes of the canopies; this method is time-consuming and laborious, making it inferior to the developed algorithms based on point clouds. The CH algorithm overestimates the actual canopy volume, but its run time is the shortest (${\overline{T}}_{\mathrm{CH}}$ = 0.277s), and therefore it can be used for real-time obstacle avoidance by agricultural robots in an orchard [58]. The VB algorithm can better reflect the 3D spatial state of the canopy. The value of V_VB falls within the mid-range for volume measurement, but it also required the longest run time (${\overline{T}}_{\mathrm{VB}}$ = 4.867s) which fluctuated greatly. Therefore, the VB algorithm could be selected to calculate canopy volume off-line, and to precisely measure relevant parameters of leaves, such as leaf area and LAD [39,59]. The AS algorithm can better express the actual shape of the actual canopy. The V_AS measurement was comparatively small, allowed more accurate reconstruction of the true shape of the canopy, and the run time of V_AS was in the mid-range among the five processing times; as a result, V_AS can be used to evaluate biomass and identify tree species. The slice-related problem is described below.

3.2.5. Discussion on the Different Number of Slices Used with CH and AS Algorithms

Different numbers of slices will produce different results for the reconstruction efforts and result in different volumes of the same C. grandis var. Longanyou canopy; furthermore, the number of slices will also affect the run time of algorithms. The slice spacing was selected as 10 and 5 cm by Fernández-Sarría and Bi, respectively [41,42]; the $\overline{H}$ of C. grandis var. Longanyou was 3.302 m combined with the above. Therefore, five different slice numbers (N = 10, 20, 30, 40, 50) were selected to discuss the influence of different numbers of slices on CH and AS algorithms from the two perspectives of calculated volume and run time.

Based on the original CH algorithm, the percentages of holes and gaps removed by different number of slices (N = 10, 20, 30, 40, 50) were 12.2%, 15.7%, 17.6%, 18.7%, 19.7% respectively. As shown in Figure 11a–e, V_CHBS decreased with an increasing N; the difference between V_CH and V_CHBS increased slightly, which resulted in the regression coefficient gradually moving away from 1, and R² decreased to a small extent. According to the above, this occurred because the method of slicing divides the total crown volume into a certain number of sub-polyhedrons along the Z-axis, removing some holes and gaps in the volume calculated by the CH algorithm. As N increased, more holes and gaps were removed so that the difference in the volume increased and the actual crown shape at each height was better approximated.

In terms of run time, the corresponding run time will increase by 1 s for each increase of ten slices. However, the volume estimated by the CHBS algorithm decreased as the increase in the number of point clouds was small, and the difference between the minimum and maximum run times for the same number of slices was no more than 0.5 s (Figure 11f). This indicates that the operation time of the CHBS algorithm was stable with a different number of point clouds.

On the basis of the original AS algorithm, the percentages of holes and gaps removed by different number of slices (N = 10, 20, 30, 40, 50) were 6.8%, 13.1%, 19.2%, 25.4%, 32.0% respectively. The relationships, trends, and causes of the differences between the V_AS and V_ASBS measurements were similar to those for V_CH and V_CHBS, and will not be repeated here (Figure 12a–e). Note that when the number of slices was small (N = 10, 20, 30), the R² between AS and ASBS algorithms was greater than that between CH and CHBS algorithms (R² of AS and ASBS was above 0.98, and the maximum R² of CH and CHBS was 0.979). The R² of AS and ASBS decreased with the increase in the slice number, and the attenuation was greater than that of CH and CHBS. In addition, the regression coefficient of AS and ASBS increased with an increase in the number of slices, and the amplification was greater than that of CH and CHBS. This shows that the AS algorithm can better remove the holes and gaps in the tree canopy than the CH algorithm, which is the same as the results presented based on the principles of these two algorithms.

In terms of run time, T_ASBS was 1 s longer than T_CHBS for the same number of slices. When the number of slices increased by ten slices, the run time of T_ASBS also increased by 1 s, and its value significantly increased. The difference between the minimum and maximum run times for the same number of slices was about 1 s (Figure 12f). This shows that the ASBS algorithm had a large range of fluctuation in run times under different numbers of point clouds.

In summary, for the canopy of C. grandis var. Longanyou, adding different slices to the CH algorithm will affect its operating efficiency. The R² values between V_CH and different V_CHBS values were 0.970–0.979, which means that the V_CH can be calculated, and the different V_CHBS also can be calculated using a linear fitting model. As the N increased in CHBS, fewer holes and gaps were removed, but the time was increased. Therefore, the CHBS algorithm is not suitable to calculate the canopy volume of C. grandis var. Longanyou. On the contrary, with the increase in N in ASBS, a greater number of holes and gaps were removed. Although the time was also increased, the R² values between AS and different V_ASBS were 0.962–0.999, and different V_ASBS values were also able to be calculated by V_AS. In addition, according to the above fitting relationship of the volumes calculated by different algorithms, the R² of V_CH and V_VB was 0.898, and the R² of V_VB and V_AS was 0.876. Therefore, V_VB can be calculated by V_CH, and V_AS can then be obtained. Combined with the above analysis of slices, it can be seen that V_ASBS can be further obtained based on the obtained V_AS. In addition, the above conclusion shows that the run time of the CH algorithm was the shortest ($\overline{H}$ = 0.277). Therefore, for the canopy of C. grandis var. Longanyou, the CH algorithm with the shortest time can be used to obtain the accuracy of the AS algorithm, and even the higher accuracy of the ASBS algorithm, which provides a better approximation of the actual shape of fruit tree canopy.

4. Conclusions

The canopy reconstruction and volume calculation of Citrus. grandis var. Longanyou were realized by the MM and five APC algorithms in this article. The specific research contents are summarized as follows:

(1)

A set of processing procedures is available for the extraction of point cloud parameters of the fruit tree canopies specifically for the ridge-tillage planting characteristics of C. grandis var. Longanyou in hilly and mountainous areas, forming a theoretical system for canopy volume analysis of this species. In addition, compared with five APC algorithms, it is difficult to select an appropriate geometric MM model, and the MM method takes a large amount of time, making it inefficient. Therefore, the APC algorithms were proven to be convenient and reliable for the calculation of canopy volume.

(2)

A suitable algorithm can be selected based on the different operational needs of C. grandis var. Longanyou: the CH algorithm is suitable for real-time obstacle avoidance and path planning for agricultural robots in orchards; LAI, LAD, and other leaf-related parameters in an offline state can be calculated by the VB algorithm; and the AS algorithm can be used for biomass estimation, tree species identification, and other operations in orchards. Furthermore, the AS algorithm was optimized by slicing in this study, and results of the ASBS algorithm were closer to the actual canopy shape. Although the running times were higher, the R² values of the pairs of CH and VB, VB and AS, and AS and ASBS were all above 0.85, indicating good linear quantitative relationships existed among different volume values. Therefore, the high-precision volume value (V_ASBS) can be calculated via the quantitative relationships using the CH algorithm, with the shortest running time.

(3)

Changing the number of slices had little influence on the volume calculated using the CH algorithm, but had a significant influence on the AS algorithm results. Furthermore, the addition of slices affected the running efficiency of the algorithm. Therefore, for the calculation of C. grandis var. Longanyou canopy volume, the CH algorithm cannot be optimized with adding slices. By comparison, the implementation of slicing optimization with the AS algorithm is dependent on specific conditions. In addition, because R² values between CH and different CHBS, AS, and different ASBS algorithms were greater than 0.95, the different volumes determined using the CHBS and ASBS algorithms under different slices can be calculated using only the CH and AS algorithms.

Discrete surfaces will be formed when α is set to a value that is too small, resulting in incomplete reconstruction results. Consequently, the effective implementation of the ASBS algorithm results in significant relationships with the distribution of the point cloud and the dynamic matching of the values of α and N. This will continue to be explored in depth in future work.