A Novel Approach for Retrieving Tree Leaf Area from Ground-Based LiDAR

Leaf area is an important plant canopy structure parameter with important ecological significance. Light detection and ranging technology (LiDAR) with the application of a terrestrial laser scanner (TLS) is an appealing method for accurately estimating leaf area; however, the actual utility of this scanner depends largely on the efficacy of point cloud data (PCD) analysis. In this paper, we present a novel method for quantifying total leaf area within each tree canopy from PCD. Firstly, the shape, normal vector distribution and structure tensor of PCD features were combined with the semi-supervised support vector machine (SVM) method to separate various tree organs, i.e., branches and leaves. In addition, the moving least squares (MLS) method was adopted to remove ghost points caused by the shaking of leaves in the wind during the scanning process. Secondly, each target tree was scanned using two patterns, i.e., one scan and three scans around the canopy, to reduce the occlusion effect. Specific layer subdivision strategies according to the acquisition ranges of the scanners were designed to separate the canopy into several layers. Thirdly, 10% of the PCD was randomly chosen as an analytic dataset (ADS). For the ADS, an innovative triangulation algorithm with an assembly threshold was designed to transform these discrete scanning points into leaf surfaces and estimate the fractions of each foliage surface covered by the laser pulses. Then, a novel ratio of the point number to leaf area in each layer was defined and combined with the total number of scanned points to retrieve the total area of the leaves in the canopy. The quantified total leaf area of each tree was validated using laborious measurements with a LAI-2200 Plant Canopy Analyser and an LI-3000C Portable Area Meter. The results showed that the individual tree leaf area was accurately reproduced using our method from three registered scans, with a relative deviation of less than 10%. Nevertheless, estimations from only one scan resulted in a deviation of >25% in the retrieved individual tree leaf area due to the occlusion effect. Indeed, this study provides a novel connection between leaf area estimates and scanning sensor configuration and supplies an interesting method for estimating leaf area based on PCD.


Introduction
Leaf area is very important for scientists in their understanding and modelling of gas-vegetation exchange processes, such as photosynthesis, evaporation and transpiration, rainfall interception, carbon flux and nutrient cycle, and is also critical to evaluate biological parameters of plant or forest.Leaf area can be estimated in situ using "direct" or "indirect" methods.
Direct methods of leaf area estimation involve manually counting and measuring leaves with a leaf area meter or image scanner [1] and using image analysis software.In some studies, the correlation between petiole length and leaf weight has been used to deduce leaf area [2].These methods are time consuming and inefficient when working with a large number of samples.Therefore, "indirect" non-contact measurements have been developed and are typically preferred.These techniques include a multitude of non-contact leaf area index (LAI) measurement tools, such as hemispherical photography, the Hemi-view Plant Canopy Analyser [3] from Delta-T Devices, the LAI-2000 Plant Canopy Analyser [4] from LI-COR Biosciences, the LP-80 AccuPAR ceptometer [5] from Decagon Devices, a mobile sensor combined with the canopy light transmittance model [6], the prototype sun-fleck LAI instrument named Tracing Radiation and Architecture of Canopies TRAC [7] and a Moderate Resolution Imaging Spectra-radiometer [8].In addition, light detection and ranging technology (LiDAR), which has been applied to airborne [9] and terrestrial instruments [10], has been used to assess canopy structure and estimate vegetation attributes [11].In recent years, LiDAR-based leaf area retrieval has become an appealing concept due to its ability to capture the structural information of canopies as 3D point cloud data (PCD).The PCD acquired using this technique records geometrical information, reflected intensity, wavelength and colour.This technique has greatly improved the accuracy of leaf area estimation and has been used in numerous studies to retrieve leaf area from various vegetation types [12].
Existing methods using terrestrial laser scanner (TLS) data to estimate leaf area can be divided into two groups: (1) Many statistical models and forestry formulas, including the Poisson gap model [13], Warren Wilson's contact frequency method [14], Monsi-Saeki theory of light attenuation following Beer's law [15], the Canopy Semi-analytic P gap and Radiative Transfer (CanSPART) model [16] and the gap size distribution method [7], have been developed to analyse the light transmission properties of canopies and reflect the spatial distribution of foliage area in tree crowns using a computer.Moreover, different space partitioning methods and projection strategies, such as point quadrat analysis [17], voxelization of 3D scattered laser points [18], and different projection methods for transforming 3D PCD into 2D hemispherical photography images [19,20], were used to retrieve tree leaf area density and canopy LAI values; (2) Computer graphics and machine vision theory were also applied for point cloud processing.Many computer techniques, such as multi-angle photographs and combined multi-view stereo vision algorithms [21,22], have been used to produce high spatial resolution 3D point clouds of vegetation structure.Some computer software, including Arbaro [23], VegeSTAR [24] and Speedtree [25], has been designed to automatically create geometric explicit models of canopy structure.Moreover, many computer algorithms for predicting canopy variables, branching structure and biomass have been proposed, including computer graphics algorithms based on morphological traits [26], neighbour-relations and local surface fitting algorithms [27], and the 3-D convex hull algorithm [28].
Although numerous studies have been performed using tree leaf area measurements based on TLS scanning data, four problematic issues arise when acquiring forestry parameters from scanned data: (1) The need to design automatic recognition algorithms for the classification of various tree organs based on the point clouds of the trees with complex branch morphology and large numbers of overlapping leaves; (2) The need to remove existing noise in PCD and build a compensation mechanism for the occlusion effect because the scanning data always suffer from a significant shadowing effect (the occlusion of canopy elements by others from the viewpoint of the LiDAR instrument) and are altered by the presence of mild wind [29,30]; (3) The need to retrieve leaf area from PCD with different spatial densities because empirical formulas show that the spatial resolution of scanning points is inversely proportional to acquisition distance [31,32]; (4) The need to transform discrete points to leaf surfaces when calculating leaf area because scanned data consist of discrete sets of points rather than surface area.
Focusing on the above questions that have plagued leaf area estimation using terrestrial LiDAR observations, the overall purpose of this study is to apply the following strategies to develop a systematic method for accurate leaf area estimation from PCD.
(1) Quantitative identification of photosynthetic and non-photosynthetic components is the primary step in effectively retrieving leaf area.Here, a series of PCD features, such as shape, orientation, structure tensor and normal vector distribution, are proposed with the semi-supervised Support Vector Machine for automatically recognizing different plant organs from plant canopy scanned data.(2) To alleviate the perspective occlusion effect and avoid the variation in spatial resolution of PCD due to varying acquisition distance, multiple scanning positions around target tree are employed to acquire full coverage of the PCD, and whole scanning data are innovatively separated into different layers to approximate the spatial resolutions of the scanned points within each layer.Moreover, the moving least squares (MLS) method is adopted to eliminate ghost points caused by blowing wind during the scanning process.(3) Following the conversion from discrete points into leaf surface data, which is aided by a novel 3D triangulation strategy with assessable thresholds, a process is conducted to automatically distinguish the gap between the leaves and partial leaf area covered by the laser pulses.
In addition, the ratio of point number to foliage area is originally determined to retrieve the leaf area for each layer.
The overview of the different steps of our method is given in the following paragraphs and schematically represented in Figure 1.

Data Collection by TLS
The experimental trees were selected on the campus of Nanjing Forestry University (32 • 08 N, 118 • 81 E) and included many well-isolated individual trees, such as michelia trees, sakura trees and Fatsia japonica.Other experimental subjects, such as rubber trees, were located in Danzhou City (19 • 30 N, 109 • 29 E), which is the largest rubber production base on Hainan Island, China.All object trees with broad leaves and canopies with varying densities were scanned using a TLS (Leica C10) and different scanning patterns.Figure 2a shows information regarding the acquisition of tree data using the scanner.Assuming that the shape of the target tree crowns is ellipsoidal, each scan was obtained using an azimuthally symmetric location with the target tree located in the centre of the experimental plot.Each TLS was placed in turn on different lateral sides and 3 m from the target tree.After the TLS scanning procedure was finished, each scan from a different angle was integrated into a single coordinate system by using a registration process to acquire full coverage scanning data from the objective trees (see Figure 2b).For Fatsia japonica and rubber tree, we used one side-lateral scan to obtain data and assess the impacts of the occlusion problem (see Figure 2c).In addition, we manually counted the number of leaves and measured leaf area using a LI-3000C Portable Leaf Area Meter.The obtained results serve as a benchmark for verifying our method results.

Reference Values Obtained Using Other Methods
To demonstrate the validity of our method, all target trees were measured using an LI-3000C portable area meter and LAI-2200 plant canopy analyser [33] for single leaf area and LAI generation, respectively.Figure 3a,b illustrates data collection by LI-3000C and image processing for retrieving leaf area.All instruments were newly procured and factory calibrated.In addition, the labour required to quantify the number of leaves in each canopy was high.For rubber tree, 30 rubber stands with areas ranging from 1.16 to 14.48 were randomly sampled [34].Eight of the 30 sampled stands were used for long-term monitoring by the Danzhou Investigation and Experimental Station of Tropical Crops, Ministry of Agriculture, Danzhou, China.The LAIs of the eight validation stands, which were measured from March 2012 to July 2014 by LAI-2200, are shown in Figure 3c.Due to the intensive labour requirements and destructive nature of this method, the lower part of one rubber tree enclosed, shown by the blue line in Figure 3d, was chosen as the experimental data.The descriptive variables for the lower part of the measurement tree were height, 10.0 m; total area of leaves in the lower part obtained by collecting deciduous leaves, 37.55 m 2 ; canopy projection area, 15.51 m 2 ; scanning point number in cloud, 166,062; tree age, 16 years; survey date of the target tree, 18 April 2014; and deduced LAI, 2.42 m 2 •m −2 .The descriptive variables for the upper part of the rubber tree were: height, 10.0 m; number of scanning points in the cloud, 104,326; and relevant LAI in proportion to the number of scanning points, 0.83 m 2 •m −2 .The whole tree was scanned from one position 2 m above the ground; thus, the scanned pattern from the bottom to top led to coarser spatial resolution data and data loss in the upper part.Thus, the LAI of the whole rubber tree determined by measuring deciduous leaves was nearly 4 m 2 •m −2 , which is consistent with the average measurement results using LAI-2000 and labelled by the red square in Figure 3c.Generally, the measurement results obtained from these methods were considered as a benchmark for evaluating the effectiveness of our proposed method based on TLS data.

Tree Organ Classification
The unknown degree of wood components likely introduces errors into leaf area estimates [35].Separating the contributions of leaf and non-leaf elements is important for our method.The detailed information acquired directly from the TLS includes colour, spatial position and intensity, but these features are influenced by the external environment and cannot reliably discriminate foliage from unripe or green branch varieties.To derive descriptions of vegetation structure at various levels of detail, some PCD features, such as shape, orientation and structure tensor, were employed for the classification of different tree organs.
For shape features, we employ a variant on the "feature saliency" for point wise analysis of PCD.Intuitively, we examine whether the local neighbourhood forms a sphere (compact and spatially isotropic), a line (a distribution with a single dominant axis of spatial variation), or a thin flat sheet (two axes of variation).These mathematical concepts are used to directly map the categories of plant structures, such as fruits (which tend to be spherical), branches (roughly linear) and leaves (largely flat) [36].The specific steps of this procedure are described as follows.
For a point p i = (x i , y i , z i ) T in the point cloud P, P ⊂ R 3 , the neighbours of point p i within the radius r N were defined as p j = x j , y j , z j T and satisfied condition ||p j − p i || ≤ r N .
The covariance matrix C P of the point p i neighbourhood is defined as where u = (1/N) ∑ N j=1 p j denotes the mean of p i neighbours.Here, we utilized eigenvalue decomposition to define a new coordinate system in which axis directions are given by eigenvectors and point variances along axes are given by corresponding eigenvalues.Let e k,i be the eigenvector of C P and λ k,i be the corresponding eigenvalue, where k = 0, 1, 2 and λ 0,i ≤ λ 1,i ≤ λ 2,i .λ k,i quantitatively shows the data variance along the axis e k,i .e 0,i = e ix , e iy , e iz corresponds to the smallest eigenvalue λ 0,i and approximates the normal vector at point p i .Next, we normalize the eigenvalues of the covariance matrix corresponding to every point p i , i.e., λ k,i = λ k,i / 2 ∑ k=0 λ k,i .Structure tensor features [37], including structure tensor planarity, structure tensor "eigenentropy" and structure tensor linearity, were achieved and are shown as Moreover, we calculated the normal vector distribution of every p i using the following equation.
where e denotes the mean normal vector of the p i neighbourhood, e 0j is the normal vector of the p i neighbour, and e = (1/N) ∑ N j=1 e 0,j .Notably, the relative magnitudes of the three eigenvalues {l 0 , l 1 , l 2 } of V p satisfying l 0 ≥ l 1 ≥ l 2 can reflect the geometrical properties of the p i neighbourhood.When l 0 ≈ l 1 ≈ l 2 reflects the isotropic spatial distributions (corresponding to fruits), l 0 ≥ l 1 ≈ l 2 indicates predominantly linear distributions such as branches and l 0 ≈ l 1 ≥ l 2 shows roughly planar distributions (leaves).
Using the above analysis, a series of features about each point p i was obtained to compose feature vectors as f p i = e ix , e iy , e iz , c 0,i , c 1,i , c 2,i , l 0,i , l 1,i , l 2,i , where e ix , e iy , e iz is the normal vector, {c 0,i , c 1,i , c 2,i } represents the structure tensor, and {l 0,i , l 1,i , l 2,i } describes the distribution of the point normal vector.
Based on the proposed point features f p , the semi-supervised Support Vector Machine (SVM) method was provided for PCD classification.Then, we manually selected 5%-10% of the points as training samples with manually labelled categories (branch or leaf).For the remainder of the unspecified PCD data, we adopted a repeated sub-sample cross-validation method based on SVM to finish the classification.

Total Leaf Area Retrieval
After using our classification approach, the PCD from the leaves were extracted from the scanning data, and the PCD from non-photosynthetic tissues were removed.In this section, we proposed using a new leaf area retrieval method combined with a novel indicator to specify the relationship between the leaf area and scanning point number.Once the indicator was defined, the leaf area for the whole canopy can be assessed.A detailed description of the method is introduced in the following sections.

Spatial Resolution of the Canopy Scanning Data
During the scanning process, the conditions of the external environment and the parameter regulation of scanning devices are uniformly preserved; however, the scanning results show that the spatial resolution of the scanned points is inversely proportional to the acquisition distances.In our experiments, the Fatsia japonica and rubber tree were represented by a single scan, which was separated into three layers according to the distance between the plant and scanner, as shown in Figure 4a.The point cloud was characterized by similar spatial resolution and laser coverage in the same layer with different properties in different layers.The achieved point cloud resolution is usually higher and the point spacing is usually smaller for small ranges, and vice versa for large ranges.If multiple stationary scanners are used in three symmetric positions around target trees, such as the michelia tree and sakura trees, the desired full coverage of the PCD was ensured and a minor occlusion effect occurred at the outer part of the canopy; however, the shadowing effect increased gradually and the point cloud resolution decreased in the centre of the canopy (see Figure 4b).According to different scanning patterns, we separated the entire canopy into several layers that suffered from different occlusion effects and had various point resolutions.Then, methods using statistical analysis to calculate the ratio of point number to foliage area in each layer were developed to deduce the total area of leaves within the canopy based on the scanning data.

Analytic Dataset Selection and Noise Removal
After performing point-wise classification using our method, the point cloud of the target trees was identified as either leaf or wood.Consequently, the leaf points isolated from the scanning data were divided into several layers according to different scanning pattern, for the scanning pattern of one scan, the canopy is divided into three layers with equal width from the nearest and furthest leaves to the scanner (Figure 4a).For the three co-registered scans around the target tree, the canopy is divided into two annulus zones of equal radius from the centre and residual edge area (Figure 4b).Thus, approximate spatial resolution is guaranteed for the scanned data in each layer.Moreover, approximately 10% of the foliage points were manually chosen as the analytic dataset (ADS) to determine the relationships between the scanning point number and foliage area.Four ADS groups were chosen from each canopy, and each group is shown in purple in Figure 4. Different ADS strategies were adopted according to various scanning patterns.Observed from the top view, two lines from the scanner and an arc enclosed one ADS of one scan, and two lines from the canopy centre and an arc enclosed one ADS of three registered scans.This strategy ensures that one ADS contains PCD with different spatial resolutions caused by variations in distance between the leaves and scanner, providing reasonable predictions for further leaf area estimates.
Due to the interference of the external environment during the scanning process, some leaves swaying in the wind are partially hit by the laser pulse, which produces some noise and ghost points located somewhere around the targets.The MLS [38] method was adopted to filter out noise and produce actual scanning data and is described below.
If the chosen ADS is A = p i ∈ R D , D = 3 , i ∈ {1, 2, ..., n}, where p i (x i , y i , z i ) is the 3D position of each point.Due to the interference or foliage jitter caused by wind, the target leaves always moved up and down during the scanning process; thus, the z i value of the scanned points varied.The goal of the MLS is to filter out deviation and find an actual z i (x i , y i ) for each point to substitute for the scanned value.Thus, the noise or deviation among the scanning data could be removed, and new fitting points p i x i , y i , z i were used to calculate the ratio in the next step.

Total Leaf Area Retrieval by Ratio ρ set l p
After removing the foliage noise points by using MLS, the Delaunay triangulation algorithm was used to transform discrete ADS points into numerous small triangles.Each triangle area was calculated by computing the determinant value of the following matrix T j = 0.5 p 1 p 2 × p 1 p 3 , where p 1 , p 2 and p 3 are the three vertexes of each spatial triangle following MLS filtering, p 1 p 2 is the vector from point p 1 to p 2 in 3-D space, and × represents a cross product.The perimeter L j of each triangle was considered as a criterion.When L j was larger than the threshold ς, the local region had a lower spatial sampling resolution than normal; thus, these triangles should be regarded as belonging to the occlusion region or the gaps between leaves.Meanwhile, the remaining triangles with perimeters less than the threshold ς represent the actual foliage surface contacted by the laser pulses, and summing their area results in the actual scanned foliage area in the ADS.The detailed schematic representations were illustrated in Section 3.2.Then, the key parameter of ratio ρ ADS layer , which indicates the ratio of the number of scanned points to the corresponding leaf area in each layer, was assessed because it is sensitive to reflecting the spot spacing versus the acquisition distance and occlusion effect.Equation ( 2) provides a description of this analysis, where num ADS layer is the number of scanned points in the ADS in each layer.
Combined with the calculated ρ ADS layer and total number of leaf points in each layer, we directly retrieved the total area of leaves L total area in the canopy by using Equation ( 3), where h is number of the separated layers in a canopy.

Determination of the Threshold ς
The threshold ς in our algorithm is very significant for the final individual tree leaf area estimate.In general, the threshold ς should "provide balance between the use of a large enough triangle perimeter to satisfy full representation of the scanned leaf surface and a small enough triangle to correctly identify the occluded areas and gaps between leaves."Here, an original approach to appraise the value of the threshold ς was proposed for all broad-leaved tree species with different resolutions and various acquisition distances.The TLS acquisitions were modelled according to the technical specifications of Leica C10, and the corresponding minimal angle resolution (increment) of the different scan resolutions is listed in Table 1.For example, under the coarse-resolution scanning conditions, the sample spacing was predefined as 0.2 m at an acquisition distance of 100 m, and the corresponding PCD number was 8 × 8 when the view field was set as 1 • × 1 • (horizontal × vertical); thus, the minimal interval angle α of the laser beams could be calculated as arctan(0.2/100)as 1 • /8, which is nearly 0.12 • Meanwhile, for every scan with a different resolution, we know the minimal angle increment τ and the distance d 1 from the scanner to the target leaves in the canopy.Thus, the Law of Sines can be used to express the corresponding sample spacing a 1 of two neighbouring laser beams projected onto the leaf surface that is perpendicular to the direction of the incident laser beams.However, the actual spacing of the sampling points varies with changes in leaf orientation and direction relative to the incident laser beam.As shown in Figure 5, θ 1 , which is between 0 • and 90 • , is the intersection angle of the normal vectors of the actual leaf surfaces and the incident laser beam.The distributions of the inclinational and azimuthal angles of the leaves in the canopy are random.We assume that the angular distribution of foliage elements is evenly distributed and set θ 1 ≈ 45 • ; thus, the sample spacing on the leaf surface was converted into the distance b 1 , which was calculated as follows: when the new sampling spacing b 1 is obtained, the average length of one triangle side has been assessed.Then, we provide a trust threshold value of ς = 3 × b 1 by performing various parameter adjustments of the scanning devices.

Plant Organ Classification
After applying our classification method to the preliminary scanned datasets of various experimental trees, we acquired descriptive properties of the PCD features and randomly chose minor portions of the points as training data for each class.Then, sub-sample cross-validation combined with the semi-supervised SVM algorithm was used to classify these scanned data, and the manually classification results were considered as reference data for validating the computer-based classification.
The detailed results are shown in Figure 6 and demonstrate the classification results of the sakura tree (Figure 6a,b), Fatsia japonica (Figure 6c) and michelia tree (Figure 6d), respectively.Classification accuracies of ≥ 90% were achieved for the point cloud of the sakura tree, but typical misclassification existed in the Fatsia japonica and michelia tree datasets, and some points belonging to thick and large branches were misclassified in regions where the local geometry was similar to the leaf surface.In contrast, parts of the leaf surfaces were classified as branches, especially at the junction between the leaves and branches or at the borders of the leaves.Quantitative assessments of the classification accuracies of each class of various target trees are recorded in Table 2.For different target trees with various organ geometries and spatial sampling, the optimal performance of the classification algorithm was affected by tuning parameter r N .Combined with the series of PCD features proposed here, our method is robust enough to handle changes in lighting and colour as the plant matures.Thus, a better foundation was established to present the applicability of leaf area appraisal and non-destructive vegetarian measurements.

Foliage Reconstruction and Leaf Area Appraisal
After classifying via our method, the classified leaf points in the canopy were used as preliminary data, among which 10% of the leaf points in each layer were randomly chosen as ADS.MLS and Delaunay triangulation were adopted to appraise the point spatial resolution versus acquisition distance and the ratio of the point number to leaf area in the cloud to deduce the total area of leaves in each layer.Specific applications of our methods for different tree species, including sakura trees and Fatsia japonica with different plant morphologies (e.g., leaf shape, leaf length, canopy architecture and growth characteristics) are demonstrated in Figure 7.Moreover, different scanning patterns, including single scan data or the registration of three scans, were used to verify the occlusion effect of our method when retrieving individual tree leaf area.
Figure 7a,b repeatedly illustrates the segmentation results of two tree species.The branch regions presented in black were accurately extracted, and the remaining regions presented in green represent foliage elements in the canopy.Meanwhile, a slicing method was developed by defining several cylindrical lateral surfaces with different radii and the same axis to divide the canopy into two or three layers, as described in Section 2.4.2.Next, MLS was used to smooth and resample some of the leaf point cloud, which was selected as ADS.3D visualization is shown in Figure 7c,e in which the initial scanning points are presented in green, and the surface-fitting operation results using the MLS method are presented in red.Because of its global convergence, the MLS method can be used to search for multiple local optimal solutions and to obtain optimal solution aggregates of the equations; thus, the MLS in this application can describe the geometry and reflect the local ecological deformation of the foliage and filter out true outlier points caused by external environmental interference.Consequently, Delaunay Triangulation was used to convert these processed point sets into leaf surface in the form of a triangular mesh.Each triangle with an edge built by connecting the nearest neighbouring points constitutes the foliage surface.The specified algorithm diagrams of the triangular mesh representation are shown in Figure 7d,f,g where the green triangles with larger perimeters outline the occluded regions of the leaves or gaps within the canopy and the blue triangles with small perimeters correspond to the local foliage surface contacted by laser pulses.The summed area of the total triangular meshes in blue reflect the local foliage surface area of the ADS exposed to laser scanning, from which the ratios of the number of points to the leaf area in each layer (ρ ADS layer ) were deduced.Combined with the total number of scanned points, the total area of leaves in each layer was easily retrieved using the methods described in Section 2.4.3.To evaluate the performance of our methods for different tree species, we randomly chose points for the ADS from scanning data of the canopy to estimate the ρ ADS layer ratio and leaf area using regression analysis.The optimal threshold values of every target tree were also deduced by using our method described in Section 2.4.4.Detailed descriptions of our results and a comparison between our methods and labouring measurements are shown in Table 3.For the Fatsia japonica and rubber tree, only a one-sided lateral scanning pattern was adopted and the canopy was partially covered by laser pulses.The results show that the ratio of ρ ADS layer of each layer decreases as the acquisition range increases (i.e., ρ ADS layer=1 > ρ ADS layer=2 > ρ ADS layer=3 ), which indicates that the occlusion effect was somewhat accentuated and that the spatial sampling decreased for the more distant portions of the canopy.For the sakura and michelia trees, a multi-location scanning pattern was adopted, the occlusion effect was effectively alleviated and evenly distributed ratios of zones at the centre and residual edge areas were achieved.For sakura tree 2, a very dense canopy makes the ρ ADS layer ratio reached a minimum value at the inner portion (layer 2) with occlusion dominating and increasing again near the edge zone (layer 1).In addition, the descriptive statistics of the total leaf area for the five individual canopy variables are also summarized in Table 3.The deviation of individual tree crown leaf area estimates between our method and the manual measurements was relatively high for the canopy data obtained from one scan and relatively low for the canopy data obtained from the multi-scanning pattern.For example, the deviation of Fatsia japonica peaked at 27.92% due to a serious occlusion effect.When three scans were combined for the complete coverage of the target trees, the approximate average deviation was 7.50% for the sakura trees.By manually counting the number of leaves in each canopy, the appraised average areas of one leaf obtained using the two methods are listed in Table 3, which indicate that our approach obtains satisfying individual tree leaf area estimates.Meanwhile, we concluded that the occlusion effect discrepancies caused by various scanning patterns has considerable impact on the overall accuracy of the results.Figure 8a depicts the measured vs. predicted individual tree leaf area values obtained using our method and manual measurements.The data points of the predicted versus measured results of different target trees followed the 1:1 line more closely, and such superior performance of our method was noted consistently in the canopies of the sakura trees 1 and 2 and the michelia tree, which were scanned at multiple locations.For the individual tree leaf area retrieval from one scan, such as Fatsia japonica and rubber tree, the predicted value was lower than the measured value because the shadow effect resulted in incomplete information acquisition.Compared with the manually measured leaf area per tree, the results estimated from our method were usually affected by the occlusion effect, ADS selection and the threshold ς value assignment, and these factors yielded wider prediction intervals for all five canopy variables.Both the ADS selection and threshold were important for the final individual tree leaf area and average area of single leaf (AASL) estimates.For each target tree, we randomly selected five ADSs with 15 different thresholds to calculate the final leaf area per tree.Combined with the known number of leaves, the variations of the deduced AASL with different thresholds are shown in Figure 8b.We found that the threshold is sensitive to the final individual tree leaf area estimate.The light blue areas represent the optimal threshold range defined by our method.
In these areas, the results showed that our obtained results for the canopies with three-scan coverage captured the variations of the manually measured leaf area well.However, for the canopy represented by only one scan, our method with optimal thresholds tended to underestimate the total leaf area of the rubber tree.This result indicated that the occlusion of objects within the canopies created a serious local deficit of material; thus, a larger threshold of deviation from the normal level obtained better results, which compensated for the loss caused by occlusion.

Comparison with Existing Methods
Our first objective was to develop an impressive and novel method for retrieving leaf area by using TLS scanning data.The core contributions of our method include the ability to address leaf area estimates by applying computer graphics, which differ from many existing methods, such as the gap probability (Pgap) model [33], the voxel-based approach [39,40], the radiative transfer model [41] and the statistical estimator of vegetation parameters [42].Specifically, Zheng et al. [40] constructed the 3D convex hull that contained the original points within one voxel, and one-half of the total surface area of this 3-D convex hull was considered to represent the actual foliage area.However, unfavourable leaf area estimate results always occurred when the points on the leaf surfaces belonging to several leaves existed in the same voxel or when some ghost points were among the scanning data (Figure 9a).Indeed, the construction of triangles to compose the partial surface of each leaf and the summation of blue triangles as the real foliage area in voxels is a superior concept (Figure 9b).According to the distance in each voxel between neighbouring laser pulses emitted with a given angular separation, Béland et al. [41] determined the number of return laser pulses trigged by leaf points in each voxel.A parametric model based on this principle yielded the leaf area density for each voxel and the total leaf area in the tree crown.Because the scanned points cannot represent the leaf surface and because the leaf point density varies in different voxels, the results of this approach based on the line-point intersection algorithm [41,43] were greatly affected by the algorithm parameters, such as the laser pulse incidence angles, distance between neighbouring laser pulses and spot spacing in the voxel (Figure 9c).After transforming discrete points to leaf surface using our triangulation processing, we were able to accurately calculate the amount of light captured by the canopy (Figure 9d).Moreover, some anomalies, including the occlusion effect among foliage elements [44], the point density variances in the canopy, the noise scattered from objects due to wind interference and the unclear representation of wood and foliage components [45], are present in the original scanning data and complicate the retrieval of leaf area data.Our research is in agreement with the strategy of [44], in which scanning data should be divided into various layers based on the acquisition distance, and occlusions mainly occur in the far layers of the crown because the scanner only provides a better view of the near layers.Moreover, similar methods for crop growth parameter estimation have been developed [21,23], of which the idea of characterizing the morphological traits of crop leaves through computer graphics theory demonstrates the promise for estimating leaf area.We further incorporated these findings and used three scans around the target canopy to realize data acquisition with full coverage.The hierarchical strategy according to the acquisition distances of the scanner was designed to process the scanned data using uniform spot spacing in each layer and by combining the strategy with the MLS method to remove the noise caused by wind.In addition, Delaunay triangulation with the predictable perimeter threshold was first proposed to transform discrete scanning points into the 3D foliage surface and calculate the ρ ADS layer ratio of the point to leaf area for deducing the total area of leaves in each layer.The results of quantifying the total leaf area in the canopy are especially relevant to manual direct measurements.

Recommendations
The following points must be considered to apply this method to broad-leaved canopy measurements.
First, because the impacts of trunks on the total scanned data and LAI estimates strongly vary during the leaf-off and leaf-on seasons [46], distinguishing between wood and foliage return signals is important for directly retrieving leaf area.The classification accuracy is particularly dependent on the feature selection, rather than the classification algorithm itself.Our method proposes using a series of PCD features that provide density and pose invariant descriptions using the properties of differential geometry.The effectiveness of these features depends on the radius r N of one point neighbourhood, which is changed by different point cloud resolutions and canopy geometry representations.Some conditions, such as the technical specifications of the scanning instrument and the required resolution of the target trees, can be controlled.The relative radius values of different tree species are always related to these conditions, from which we can easily obtain the desired r N .Moreover, a tree is invariably composed of three parts, i.e., trunk (big branches), branch and foliage.Trunks are always much thicker than branches, which likely makes the local surface of the trunk more indicative of planar properties.Thus, we recommend dividing each tree PCD into two segments according to the height; the lower tree segment contains primarily the trunk and foliages, whereas the upper tree segment contains branches and foliages.Then, a Gaussian mixture model (GMM), SVM or Libsvm [47] with suggested parameters can be performed on each segment.Because the number of PCD features is much smaller than the number of samples, the linear kernel and Gaussian radial basis (RBF) kernel are preferred for SVM and Libsvm, and SVM is more time consuming than the Libsvm approach.Matlab is preferred because it has GMM and SVM packages for direct operation.The overall accuracy of these methods with our constructed PCD features exceeds 85%, which is near the preliminary experiments shown in [36].Moreover, the multispectral or hyperspectral information collected by spectral scanning instruments such as the RIEGL terrestrial scanner (RIEGL USA Inc., Orlando, FL, USA) could be used to improve the final classification.
Second, the spatial resolution of the point clouds decreases as the acquisition distance increases [48]; thus, it is very important to estimate the accuracy ratio ρ ADS layer , which depends on two factors, i.e., the chosen ADS and the magnitude of the threshold ς.Here, we suggested that ADS should contain at least 10%-20% of the total scanned data when using random selection.Leaf PCD at different heights and with different densities and positions should be included in one ADS.Repeated random sub-samples of ADS with cross-validation are also beneficial for predicting a reasonable ρ ADS layer distribution in each layer.Table 3 illustrates that our method yields narrow prediction intervals of ρ ADS layer for all five canopy variables.If the target tree has an irregular shape or architecture, the division of the canopy into different segments according to the distance from the scanner still can be performed.Choosing three different ADSs with cross-validation or enlarging the coverage of one ADS can be more indicative of the true magnitudes of ρ ADS layer .In the same way, the threshold ς is equally significant for the resulting estimations.In Section 2.4.4 and Figure 8b, we have stated a mean for predicting the threshold interval and shown the effect of threshold selection on the final individual tree leaf area estimate.Here, we also recommend incorporating additional thoughts, such as the ray tracing algorithm within the voxel [41] or the zenith and azimuth angle distribution of leaves in the canopy [49] for further assessing the threshold.In addition, 3D Delaunay triangulation has always seemed complex and requires a large amount of computation.The 2D Delaunay triangulation algorithm could be adopted in our work to simplify the calculation process.First, a projection strategy combined with displacement transform was designed, which assigns a displacement vector for all points in one ADS and projects the points from the spatial state p i (x i , y i , z i ) onto the plane as p i (x i + n • (z i − z min ) , y i ) or p i (x i , y i + n • (z i − z min )).Here, n represents any arbitrary integer and z min is the minimum value of the corresponding ADS.Then, the morphology structures of the projected leaf points at similar heights remain unchanged, and the overlap and congestion between leaves at different heights are removed.Second, the 2D Delaunay triangulation algorithm was created for triangle meshes based on the projected points p i , and the point indices defining each triangle that make up the triangulation were obtained.Finally, according to the point indices, a reverse projection was operated to realize 3D Delaunay triangulation for the scanned data in ADS.Performing such an approach could obviously improve the efficiency of the algorithm and possibly make 3D triangulation unnecessary for our algorithm.
Third, many factors affect the performance of our method, namely object occlusion, noise point caused by wind, optimal selection of algorithm parameter and differentiation of photosynthetic and non-photosynthetic components.The configuration of our computer is a 5th Generation Intel Core i7 processor, 8 GB RAM and the Win 7 Ultimate operating system.Both the coarser-and coarse-resolution scans are preferred when using our method due to the enormous amounts of point cloud data, resulting in a heavy computational load.However, our algorithm was able to obtain promising leaf area retrieval results for two sakura trees and one michelia tree using three scans around the target canopy.Using this scanning pattern, the deviation between the calculated and labour measured results was less than 10%.When one scan of the target Fatsia japonica and rubber trees was considered, the deviation was greater than 30% due to self-occluded areas.Although our method is not suitable for coniferous trees because needles are difficult to represent by triangle meshes, the phenotypic traits of broad leaves with wide surfaces are suitable for our algorithms.Different species have varying leaf inclination angle and leaf density distribution, resulting in different degrees of occlusion and threshold, although our method can still provide a considerable and reasonable approach for assessing the real leaf area covered by laser pulses.In addition, our method could be used to estimate the total leaf area of heterogeneous forests at the forest plot level without restricting the number of trees.Two general scanning patterns for forest measurements, i.e., terrestrial scanner alongside the forest or drone with scanner hovering over the forest plot, are commonly adopted, which are shown in Figure 10.The scanning range or the distance from the scanner to each canopy centre can be estimated through the automatic algorithm for every tree crown delineation [50].Moreover, the minimal angle increment corresponding to the different scan resolution is the basic parameter that can be easily found from the technical specification brochure for various brands and types of laser scanners.Then, the threshold ς for each tree can be extrapolated.Following the conversion of discrete laser points into leaf surface information, the ratio ρ can be calculated.Thus, our proposed method yet provides an effective total leaf area estimation for forest plots from scanned data.Each colour and circle represents various tree models resulting from tree crown delineation algorithms; therefore, the distance from the scanner to each tree can be calculated and combined with the minimal angle increment of various scanning patterns, rendering leaf area retrieval possible.(a) Top view of the experimental setup, where the scanner is alongside the forest plot; and (b) side view of the experimental setup, using drone loading with 3D scanner to collect forest PCD.

Conclusions
The novel method proposed in this study produced the true leaf area based on terrestrial laser scanner (TLS) scanning data.First, some point cloud data (PCD) features, including shape, normal vector distribution and structure tensor, can capture different spatial distribution patterns associated with different tree organs, i.e., leaves and branches.The accuracy of our method for final classification reached up to 90%.For the classified leaf points, we used the moving least squares (MLS) method to remove the noise caused by wind and the resample data for the leaf surface.Second, two scanning patterns including only one scan and three scans around the target canopy were adopted, and scanned data were separated into various layers according to the acquisition distance from the scanners.Specific layer subdivision strategies corresponding to various scanning patterns were used to determine the proximate spatial resolutions of the scanned points in each layer.Then, we used PCD triangulation to estimate the fraction of missing plots in each layer and predicted the point number to leaf area ratio from the analytic dataset (ADS).The best threshold for PCD triangulation was deduced from specific scanning parameters.The calculated ratio combined with the scanning point number in each layer was used to quantify the total true leaf area in the canopy.Finally, a large amount of work was performed using different devices to validate our results.Comparison results show that our method obtained a better recovery of 93%-95% for individual tree leaf area based on three scans around the target tree, and perspective occlusion effect decreased the recovery of individual tree leaf area to 70%-84% when the target tree was scanned from only a single position.Some issues, including object occlusion, the noise caused by wind, and miscalculation due to the non-photosynthetic part of the plants can also be considered in our work.Briefly, our method provides a novel method for appraising the true leaf area in forest stands, and has promising potential for sustainable forest management practices and future precise agricultural methodologies.

Figure 1 .
Figure 1.Schematic representation to illustrate the main steps of our method.

Figure 2 .
Figure 2. Experimental setup and scanning data collection: (a) Using terrestrial laser scanner (TLS) for various scenarios and for scanning different tree species; (b) Experimental set up in a triangular or rectangular plot with scanners placed at each side and lateral locations to ensure full coverage of the target trees, which included one michelia tree and two sakura trees; and (c) Only one scan was arranged for Fatsia japonica and rubber tree to analyse the occlusion effect existing in the LiDAR data.

Figure 3 .
Figure 3.Comparison of manual measurement means to verify our method: (a) Measured target leaf area using the LI-3000C.(b) Measured deciduous leaf area using image scanner and image analysis software.(c) Estimated variations of the mean leaf area indexes (LAIs) of the chosen rubber trees using the LAI-2200.Different colours represent the estimated LAIs of the target rubber trees with different ages.The error bars show the standard deviation.(d) The lower part of one rubber tree was selected for this study and is shown in the blue box.

Figure 4 .
Figure 4. Point clouds acquired using two scanning patterns: (a) Top view of an individual Fatsia japonica tree.The canopy was represented by one scan, which was separated into several layers according to the distance to the scanner.Each purple region contains one analytic dataset (ADS) that provides a prediction of point cloud data (PCD) density variations with changing acquisition distance.(b) Top view of three co-registered scans collected from different scanning locations and the canopies of the target trees, including the michelia, sakura and soapberry trees, were separated into two annulus zones (layers) according to the distance from the scanners to the canopy centre.Each purple region contains one ADS and includes both internal and external canopy data.

Figure 5 .
Figure 5. Schematic illustrating the relationship between the minimum angle interval and the corresponding sampling spaces.The laser beams are emitted from a point light source with an interval angle of τ: (a) the leaf surface perpendicular to the direction of the incident laser beam with sampling space a 1 ; and (b) the actual leaf surface with a random inclinational angle in the canopy.θ 1 is the intersection angle between the normal vector of the leaf surface and incident laser beam, which is assigned as the mean value around 45 • .Thus, the average sample spacing of the scanned data is adjusted to b 1 .

Figure 6 .
Figure 6.Scatter plots of point cloud data using our proposed method to automatically classify into leaves (green) and branches (brown).The obtained foliage points are prepared for leaf area estimates: (a,b) the upper and lower parts of sakura tree 1; (c) Fatsia japonica; and (d) michelia tree.

Figure 7 .
Figure 7. (a-g) Schematic diagrams illustrating the scanning data separated into several layers according to acquisition distance and leaf surface and reconstructed by using Delaunay triangulation.(a) Three scan registrations of a sakura tree separated into two layers and zones for the centre and residual edge areas.(b) One scan of a Fatsia japonica separated into three layers according to acquisition distance.(c,e) The partial foliage points (green) chosen from the scanning data of different trees as the ADS and red points represent the resampling leaf points through MLS filtering.(d,f,g)The ADSs of the sakura tree and Fatsia japonica, and triangle meshes were built based on the red points, where the green triangles with larger perimeters represent the gaps within the canopy or the partially occluded regions and the blue triangles with small perimeters represent the local leaf surface contacted by laser pulses.

9Figure 8 .
Figure 8. Comparisons of the estimated leaf areas between our algorithm and the manual method for five target trees: (a) scatter plots of the reference values obtained using LI-3000C versus the LiDAR-predicted individual tree leaf area obtained using our method; and (b) the obtained average area of single leaf (AASL) variations of different tree species with changes in threshold and different ADSs.The light blue areas labelled as optimal forecast threshold ranges defined by our algorithm for AASL prediction.

Figure 9 .
Figure 9. Schematic diagrams illustrating the superiority of our algorithm: (a) Wire-frame representation of the 3-D convex hull that contains partial scanned points belonging to several broad leaves.Obviously, one-half of the total surface area of this hull is not equal to the actual foliage area.(b) Based on the scanned points, our method can clearly delineate the foliage surface (blue triangles) covered by the laser pulses.The summation of these triangle areas allows for precise estimates of the actual foliage area in the voxel.(c) Several methods derive LAI by using a line-point intersection algorithm.However, this method is arduous for computing intersections because the distance between the neighbouring laser pulses is always inconsistent with varied point spacing.(d) After transforming the discrete scanned points into leaf surfaces, the line-surface intersection algorithm can be used to precisely determine the number of return pulses trigged inside the voxel by the leaves.

Figure 10 .
Figure10.Schematic diagrams illustrating the method applied to the forest plot measurements.Each colour and circle represents various tree models resulting from tree crown delineation algorithms; therefore, the distance from the scanner to each tree can be calculated and combined with the minimal angle increment of various scanning patterns, rendering leaf area retrieval possible.(a) Top view of the experimental setup, where the scanner is alongside the forest plot; and (b) side view of the experimental setup, using drone loading with 3D scanner to collect forest PCD.

Table 1 .
Leica C10 main technical specifications and the qualitative results of the minimum angle resolution at 100 m with different scanning resolutions.

Table 2 .
Overall accuracy assessment for the classification of different tree species.

Table 3 .
Statistics of leaf area estimates for the five individual trees using our method and LI-3000C.