# Study of Subtropical Forestry Index Retrieval Using Terrestrial Laser Scanning and Hemispherical Photography

Authors: Ting Yun; Weizheng Li; Yuan Sun; Lianfeng Xue
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

In order to retrieve gap fraction, leaf inclination angle, and leaf area index (LAI) of subtropical forestry canopy, here we acquired forestry detailed information by means of hemispherical photography, terrestrial laser scanning, and LAI-2200 plant canopy analyzer. Meanwhile, we presented a series of image processing and computer graphics algorithms that include image and point cloud data (PCD) segmentation methods for branch and leaf classification and PCD features, such as normal vector, tangent plane extraction, and hemispherical projection method for PCD coordinate transformation. In addition, various forestry mathematical models were proposed to deduce forestry canopy indexes based on the radiation transfer model of Beer-Lambert law. Through the comparison of the experimental results on many plot samples, the terrestrial laser scanner- (TLS-) based index estimation method obtains results similar to digital hemispherical photograph (HP) and LAI-2200 plant canopy analyzer taken of the same stands and used for validation. It indicates that the TLS-based algorithm is able to capture the variability in LAI of forest stands with a range of densities, and there is a high chance to enhance TLS as a calibration tool for other devices.


---

## Body

## 1. Introduction

The forest has an irreplaceable status in improving the Earth’s environment of human habitation and slowing down global environmental degradation trend. Tree canopy performs three important functions such as tree’s photosynthesis, interchange of gases and evaporation of water between the atmosphere and the plant body, and plant growth and development, so how to solve and explore the fine measurement of tree canopy and find an effective way to calculate the canopy indexes, such as leaf area index, crown density, and leaf area density, is a main task in forestry research and an important topic of inverse analysis aimed at plant growth mechanism.In recent years, there have been several methods of obtaining canopy morphology structure and measuring trees parameters. Firstly, plant canopy analyzer and traditional mechanical instruments for trees measurement are inefficient and are affected by leaf overlapping and aggregation effects in addition to the solar zenith angle. Secondly, hyperspectral radiometer has been used to measure the canopy from top to down and usually includes recorded reflection from the ground [1, 2]. It is capable of producing one-shot topographic and spectral intensity information, which will enable a simultaneous study of structural and biochemical vegetation parameters, but this method is constrained by atmospheric conditions and accuracy of aerial photography and it is not theoretically supported by modern physical models. Additionally, airborne LIDAR can be used to monitor plant biomass and growth [3, 4], but merely relying on finitely echo waveform with attenuation and noise to account for detailed forest features is far from enough. Thus, airborne LIDAR data are always collected for generating digital forestry models at national level or large scale instead of fine scale with single tree measurement [5, 6]. Moreover, the TLS scans pulse laser over the full upper hemisphere and part of the lower hemisphere by using a mirror rotating in the vertical plane (the zenith scan) and rotation of the instrument to provide azimuthal coverage. The time-dependent intensity of reflected light from each laser pulse is recorded, providing a waveform that may include responses from multiple targets [7, 8]. The laser produces an energy pulse that is distributed in time. The shape of the outgoing pulse is consistent and has a well-defined peak. Reflections from targets are time-delayed copies of this pulse where the range to the target is simply inferred from the arrival time of the pulse peak relative to the peak of the outgoing pulse. TLS can mosaic the multiview scanning data, thereby ensuring scanning data completeness and reducing the impact of foliage cover and eventually getting 3D point cloud model of real stumpage. Due to high-precision and high-density characteristics, the TLS-based method is taken as the most appropriate method for tree measurement, and it can also serve as a calibration tool for airborne laser scanning and other measurement applications with ground sampling [9, 10].In recent years, there have been primarily several methods to calculate tree indexes from TLS PCD.(1)
Different space partitioning methods and projection strategies on PCD were used to calculate vegetation architecture and foliage assemblage in each space cell and then to indicate trees’ growth index [11, 12]. For example, Zheng and Moskal [13] presented a new voxel-based method with line quadrat direction to retrieve the biophysical characteristics of the forest canopy including extinction coefficient, gap fraction, and overlapping effect along the direction of the line quadrat and estimated effective leaf area (ELA) from TLS point cloud data. Béland et al. [14] investigated the use of a voxel-based approach to retrieve leaf area distribution of individual trees from PCD and provided vertical as well as radial distributions of leaf area in individual trees to estimate savanna vegetation structural parameters. Zheng and Moskal [15] developed circular point cloud slicing to explore the spatial variation of point density for both azimuthal angular and radial directions; his result showed that comprehensive scan combination could fully represent the canopy structure and structural variation of the heterogeneous forest stand. Cifuentes et al. [16] processed the collected 3D point clouds of forest canopy using fully representative voxel-based and multiscale models. Moreover, ray tracing algorithms are designed to transform these discrete PCD into hemispherical view, in order to assess their impact on the gap fraction estimates derived from TLS data.(2)
The physical model and forestry formula were applied to PCD for forestry index calculation. For example, the semi-analytic Pgap and Radiative Transfer (CanSPART) model was designed to predict tree gap probability (Pgap) profiles. Then, a canopy structure model based on simple geometric forms and parameterized with plot-scale statistical biometric data was used to predict gap probability (Pgap) profiles. The experimental result showed that CanSPART mathematical model performed better in clumped canopies than a simple exponential model [17]. Another method combining beer laws and using least square fitting techniques to reconstruct the normal vectors is presented, so it can indirectly and nondestructively retrieve foliage elements’ orientation and distribution from PCD obtained by using a terrestrial laser scanning (TLS) approach [18]. Moreover, Jensen’s inequality and Beer-Lambert law are also adopted to construct a parametric model for estimating leaf area density at the voxel scale from TLS data. Compared with existing ray-tracing algorithms, the model, involving computational geometry methods, can retrieve more useful information on forest structure and leaf area [19].(3)
Computer graphics and computer vision theory were studied for TLS point cloud processing. To process incomplete TLS data sets, some methods based on a structure-aware global optimization approach (SAGO) or allometric theory are proposed [20, 21]. These methods obtain the approximate tree skeleton from a distance minimum spanning tree and define the stretching directions of the branches on the tree skeleton. Then, missing branch data are recovered using different similarity functions and tree skeleton is smoothed by employing Laplacian function. Finally, the additional branches are synthesized and leaves are added to form the plausible crown geometry. Other researches [22] used a series of cylinders to build up visible branches, and the major stems in a relatively dense managed forest could be located by proposed automatic stem-mapping algorithm. For the nonvisible branches, crown feature points are extracted and hierarchical particle flow technique is designed to synthesize high order branches; thus visually convincing tree models that are consistent with scan data are produced and more pieces of information of scanned leaves can be used for reconstruction [23].Although many researchers have done a lot of valuable work in agriculture or forestry measurement, how to design an effective and convenient forestry monitoring method to estimate biomass data of woods is still a research hotspot. Hemispherical photography (HP), also known as fisheye or canopy photography, is a technique to estimate solar radiation and characterize plant canopy geometry using photographs taken looking upward through an extreme wide-angle lens [24]. HP entails five steps: photograph acquisition, digitization, registration, classification, and index calculation. Registration, classification, and calculation are accomplished using dedicated hemispherical photography analysis software. According to the principal of hemispherical photography, here we designed a new method to calculate the forestry index from PCD, which also include similar steps: PCD classification to identify leaf part and branch part, using hemispherical projection to transform leaf PCD into HP pattern, and calculation forestry index based on the PCD in HP pattern. Then, in order to verify the effectiveness of our new method, we took detailed comparison in many aspects to illustrate strengths and weaknesses of each method. The rest of this paper is organized as follows. In the section of materials and methods, we firstly used full-frame fisheye camera to capture bottom views of forestry canopy and designed image segmentation methods to retrieve the gap fraction and effective LAI of the canopy. Secondly, we used Leica C10 TLS to scan the same forestry to obtain different density PCD and meanwhile presented PCD feature extraction and classification methods to figure out normal vectors of PCD and leaf inclination angles. Besides, LAI was deduced from proposed hemisphere projection and application of forestry mathematical models. In the section of experiments and results, we compared three methods in many aspects to illustrate the effectiveness of our TLS-based method and provide concluding remarks for future work.


## 2. Data Collection and Analyses

### 2.1. Using Hemispherical Photography for Forestry Analysis

We used Canon 5D MARK III full-frame digital single lens reflex and Canon EF 24–105 mm f/4 L fisheye USM Lens as the tools for capturing high-precision bottom view of forestry canopy. Figure1 shows the process of experimental data collection by our fisheye camera, and the captured hemispherical photographs of forestry stands are shown in Figure 2.Figure 1
Canopy image collection by using fisheye camera.
(a)
(b)Figure 2
Hemispherical photographs of forest canopy.
(a)
(b)We took hemispherical photographs at 10 a.m., so the photos are affected by sunlight with nonuniform brightness. This paper presents an image segmentation algorithm with strategy of adaptive threshold selection to classify pixels of leaf element from HP denoted byIfisheye in order to eliminate errors and occlusion caused by branches and background. Specific formulas are as follows:(1)Bi,j=0Bi,j≥avei,j∗eavei,j-β1α1sky0Bi,j<avei,j∗eavei,j-β2α2branch,(2)avei,j=∑i=1m∑j=1mSi,jm∗m,where Ifisheye is divided into many small subblocks B with size of m∗m,α1,β1,α2,β2 are regulation factors, and S represents saturation component of each subblock B. Through adjustment of avei,j value, the optimal image segmentation results of leaf element extraction from HP are achieved. Consequently, we divide hemispherical photograph Ifisheye into n concentric rings with the same width, as shown in Figure 2, by red lines, where n=9. Each concentric ring corresponds to the different solar zenith angle θ. Then we can calculate canopy gap fraction Tθ of each solar zenith angle using the following formula:(3)Tθ=sumBi,j,θ==0sumBi,j,θ~=111.


### 2.2. Forestry Analysis Based on PCD

#### 2.2.1. Data Collection by TLS

We used Leica Scan Station C10 to scan the same woods and obtain their PCD for processing. The specifications of C10 are described in Table1, and Figure 3 shows our experiment about tree data acquisition using Leica C10 scanner.Table 1
The specifications of Leica Scan Station C10.
C10 technical details
Accuracy of single measurement[position/distance/angular (horizontal/vertical)]
6 mm/4 mm/60μrad/60 μrad (12′′/12′′)
Range
300 m @ 90%; 134 m @ 18% albedo (minimum range 0.1 m)
Scan resolution
Spot size
From 0 to 50 m: 4.5 mm (FWHH-based);7 mm (Gaussian-based)
Point spacing
Fully selectable horizontal and vertical range; <1 mm minimum spacing
Laser plummet
Laser class: (IEC 60825-1)Centering accuracy: 1.5 mmLaser dotted diameter: 2.5 mm
Laser color
Green, wavelength = 532 nm visible
Environmental lighting
Fully operational between bright sunlight and complete darknessFigure 3
Using TLS to scan trees on our campus.
(a)
(b)The scanning PCD of woods are shown in Figure4. In order to calculate the effective leaf area index, we firstly design feature extraction methods to extract leaf PCD from the whole forestry PCD.Figure 4
The forestry PCD obtained by TLS.


#### 2.2.2. PCD Features Extraction

Reliable feature extraction from 3D point cloud data is an important problem in many application domains, such as reverse engineering, object recognition, industrial inspection, and autonomous navigation. Many researchers have proposed many kinds of algorithms to extract the geometric features from 3D point cloud data [20, 25]. Tree’s PCD obtained by laser scanner have properties of both color and space location, but under the interference of external environment, such as illumination change and leaf swaying in the wind, the existing features of PCD, such as leaf, branch, and fruit, are inadequate to identify different tree organs. Therefore, in this paper new PCD features of topological structure and tangent space information are designed for PCD classification.


#### 2.2.3. Topological Structure of PCD

We denote PCD set byP, P⊂R3, and pi=xi,yi,ziT is the point that belonged to set P. The k nearest points of pi are denoted by pj′=xj,yj,zjT with mean value u=1/k∑j=1kpj, and the covariance matrix of points pi and its neighbor is defined as Cp=1/N∑j=1Npj-upj-uT. After taking eigenvalue decomposition, such as Cpivi=λivi, the eigenvalue λi=λ0,i,λ1,i,λ2,i and eigenvector vi=v0,i,v1,i,v2,i of Cp are calculated, where the eigenvalue λi quantitatively shows the data variance along the axis vi or the compactness of the point distribution along the axis vi, so the v0,i corresponding to the smallest eigenvalue λ0,i is the normal vector of pi. Having got the norm vector of each PCD, we can use the following formula to calculate the covariance matrix Vpi of norm vector: (4)Vpi=1N∑j=1Nv0,j-v-iv0,j-v-iT,where v-i=1/N∑j=1Nvj is the average norm vector of pi’s neighbors. Consequently, we perform singular value decomposition on Vpi to identify the principal components of P’s spatial distribution. Specifically, we analyze the relative magnitudes of the top three eigenvalues ei=e0,ie1,ie2,i of Vpi. We expect e0,i≈e1,i≈e2,i for isotropic spatial distributions (corresponding to fruits), e0,i≥e1,i≈e2,i for predominantly linear distributions (branches), and e0,i≈e1,i≥e2,i for roughly planar distributions (leaves). In brief, the topological structure of PCD can be represented by ei=e0,ie1,ie2,i.


#### 2.2.4. Tangent Space Character of PCD

Next, we design algorithm to calculate local tangent space character of PCD and determinek nearest neighbors pj′ of pi, pj′=p1′,p2′,…,pk′. We assume that a set of data points are sampled from a d-dimensional affine subspace; that is,(5)pj′=ci+Qiθj+εj1≤j≤k,where εj∈R3 represents noise vector, θj∈Rd is projection coordinates about pj′ on the local tangent space, and Θi=θ1,θ2,…,θk. d≤3. ci∈R3 is the origin coordinates of the tangent space and Qi∈R3×d is a matrix which forms an orthonormal basis of the affine subspace. The problem of linear manifold learning amounts to seek ci,Qi,θj to minimize the reconstruction error E; that is,(6)min⁡ci,Qi,θj∑j=1kpj-ci-Qiθj22=min⁡ci,Qi,ΘP-ci-QΘ22.The matrix of pi neighborhood is denoted by Xi=p1′,p2′,…,pk′, and we extract local information by calculating the eigenvectors and eigenvalue of the correlation matrix Xi-p-i′eTTXi-p-i′eT. Here p-i′=1/k∑j=1kpj, and e is a k-dimensional column vector of all ones; that is, (7)XiI-1k11TXiT=UiΛiUiT,where Ui  =ui1,ui2,…,uik is orthogonal matrix, and the diagonal elements of the diagonal matrix Λi are monotone decreasing, so the local tangent space information for the sample point pi is calculated:(8)ci=1kXi1,Qi=ui1,ui2,…,uid,Θi=QiTXiI-1k11T.From the above derivation, we can calculate the column vectoruim which is corresponding to the smallest diagonal element of Λi, and uim is also the normal vectors of local tangent space on pi.After the above analysis, a series of features about each pointpi are obtained, described as Fpi=v0,i,e0,i,e1,i,e2,i,uim, where v0,i=vi,x′,vi,y′,vi,z′ is the normal vector of pi, e0,i,e1,i,e2,i represent the distribution of v0,i, and uim are normal vectors of local tangent space on pi.


### 2.3. PCD Classification Using GMM and EM Algorithm

The scanning PCD data belong to the two semantic classes (branch and leaf). In this section, we combine Gaussian mixture model (GMM) classifier and expectation maximization (EM) algorithm to classify forestry PCD.GMM is a weighted sum ofA component Gaussian densities as given by (9),(9)pFp∣λ=∑i=1AωigFp∣μi,σi,where Fp is a continuously valued data vector of PCD (i.e., measurement or features), ωi,i=1,…,A, are the mixture weights, and gFp∣μi,σi, i=1,…,A, are the component Gaussian densities. PCD can be divided into two categories, which are branch and leaf, respectively, so A=2. Each component density is a form of Gaussian function, with mean vector μi and covariance matrix σi, and ωi is the weight coefficient of each class. The expansion formula of g is defined as follows: (10)gFp∣μi,σi=12πσi1/2·exp⁡-12Fp-μi′σi-1Fp-μi.The mixture weights satisfy the constraint that∑i=1Mωi=1. The complete Gaussian mixture model is parameterized by the mean vectors, covariance matrices, and mixture weights from all component densities. These parameters are collectively represented by the notation, λ=ωi,μi,σi,i=1,…,A. Then the expectation maximization (EM) algorithm is proposed to maximize the likelihood pFp∣λ of the data Fp drawn from an unknown distribution. In each iterative process of the algorithm, two steps of EM are executed. Firstly, estimate the distribution of the hidden variable according to the current value of the parameters; secondly, modify the parameters in order to maximize the joint distribution of the data and the hidden variable. Specific formula is expressed as follows:(11)λ∗=argmax⁡λ∏j=1N′∑k=1AωkgFp∣λ,where N′ represents total number of the calculated PCD. The output of the GMM and EM generates an independent classification for each of the 3D PCD based solely upon PCD’s feature vectors and confirmed in our experiments; such classification can be quite effective and the different plant organs are distinguished; thus, the whole PCD can be divided into two types: branches Pbranch and leaves Pleaf.


### 2.4. Spherical Projection of Classified PCD

In this section, we convert the segmented woods PCD from Cartesian coordinate system to spherical coordinate system with a radius of one, which projects the PCD onto the surface of a hemisphere. The aim of the stereographic projection was to project all points from the surface of the upper hemisphere (3D space) into the projection plane that is a 2D space. For example, the PCDPleafxi,yi,zi projected onto the spherical surface would be Pleaf′ri,φi,θi, where φ is azimuth, and θ is zenith angle. The projection process is computed and converted using (12)ri=xi2+yi2+zi2,φi=arctanyixi,θi=arccosziri.After the projection transform, each pointPleaf′ri,φi,θi is projected onto the X-Y plane and the hemispherical photograph Ifisheye′ transformed from PCD is obtained. In the following section, we use forestry models to calculate canopy indicators.


## 2.1. Using Hemispherical Photography for Forestry Analysis

We used Canon 5D MARK III full-frame digital single lens reflex and Canon EF 24–105 mm f/4 L fisheye USM Lens as the tools for capturing high-precision bottom view of forestry canopy. Figure1 shows the process of experimental data collection by our fisheye camera, and the captured hemispherical photographs of forestry stands are shown in Figure 2.Figure 1
Canopy image collection by using fisheye camera.
(a)
(b)Figure 2
Hemispherical photographs of forest canopy.
(a)
(b)We took hemispherical photographs at 10 a.m., so the photos are affected by sunlight with nonuniform brightness. This paper presents an image segmentation algorithm with strategy of adaptive threshold selection to classify pixels of leaf element from HP denoted byIfisheye in order to eliminate errors and occlusion caused by branches and background. Specific formulas are as follows:(1)Bi,j=0Bi,j≥avei,j∗eavei,j-β1α1sky0Bi,j<avei,j∗eavei,j-β2α2branch,(2)avei,j=∑i=1m∑j=1mSi,jm∗m,where Ifisheye is divided into many small subblocks B with size of m∗m,α1,β1,α2,β2 are regulation factors, and S represents saturation component of each subblock B. Through adjustment of avei,j value, the optimal image segmentation results of leaf element extraction from HP are achieved. Consequently, we divide hemispherical photograph Ifisheye into n concentric rings with the same width, as shown in Figure 2, by red lines, where n=9. Each concentric ring corresponds to the different solar zenith angle θ. Then we can calculate canopy gap fraction Tθ of each solar zenith angle using the following formula:(3)Tθ=sumBi,j,θ==0sumBi,j,θ~=111.


## 2.2. Forestry Analysis Based on PCD

### 2.2.1. Data Collection by TLS

We used Leica Scan Station C10 to scan the same woods and obtain their PCD for processing. The specifications of C10 are described in Table1, and Figure 3 shows our experiment about tree data acquisition using Leica C10 scanner.Table 1
The specifications of Leica Scan Station C10.
C10 technical details
Accuracy of single measurement[position/distance/angular (horizontal/vertical)]
6 mm/4 mm/60μrad/60 μrad (12′′/12′′)
Range
300 m @ 90%; 134 m @ 18% albedo (minimum range 0.1 m)
Scan resolution
Spot size
From 0 to 50 m: 4.5 mm (FWHH-based);7 mm (Gaussian-based)
Point spacing
Fully selectable horizontal and vertical range; <1 mm minimum spacing
Laser plummet
Laser class: (IEC 60825-1)Centering accuracy: 1.5 mmLaser dotted diameter: 2.5 mm
Laser color
Green, wavelength = 532 nm visible
Environmental lighting
Fully operational between bright sunlight and complete darknessFigure 3
Using TLS to scan trees on our campus.
(a)
(b)The scanning PCD of woods are shown in Figure4. In order to calculate the effective leaf area index, we firstly design feature extraction methods to extract leaf PCD from the whole forestry PCD.Figure 4
The forestry PCD obtained by TLS.


### 2.2.2. PCD Features Extraction

Reliable feature extraction from 3D point cloud data is an important problem in many application domains, such as reverse engineering, object recognition, industrial inspection, and autonomous navigation. Many researchers have proposed many kinds of algorithms to extract the geometric features from 3D point cloud data [20, 25]. Tree’s PCD obtained by laser scanner have properties of both color and space location, but under the interference of external environment, such as illumination change and leaf swaying in the wind, the existing features of PCD, such as leaf, branch, and fruit, are inadequate to identify different tree organs. Therefore, in this paper new PCD features of topological structure and tangent space information are designed for PCD classification.


### 2.2.3. Topological Structure of PCD

We denote PCD set byP, P⊂R3, and pi=xi,yi,ziT is the point that belonged to set P. The k nearest points of pi are denoted by pj′=xj,yj,zjT with mean value u=1/k∑j=1kpj, and the covariance matrix of points pi and its neighbor is defined as Cp=1/N∑j=1Npj-upj-uT. After taking eigenvalue decomposition, such as Cpivi=λivi, the eigenvalue λi=λ0,i,λ1,i,λ2,i and eigenvector vi=v0,i,v1,i,v2,i of Cp are calculated, where the eigenvalue λi quantitatively shows the data variance along the axis vi or the compactness of the point distribution along the axis vi, so the v0,i corresponding to the smallest eigenvalue λ0,i is the normal vector of pi. Having got the norm vector of each PCD, we can use the following formula to calculate the covariance matrix Vpi of norm vector: (4)Vpi=1N∑j=1Nv0,j-v-iv0,j-v-iT,where v-i=1/N∑j=1Nvj is the average norm vector of pi’s neighbors. Consequently, we perform singular value decomposition on Vpi to identify the principal components of P’s spatial distribution. Specifically, we analyze the relative magnitudes of the top three eigenvalues ei=e0,ie1,ie2,i of Vpi. We expect e0,i≈e1,i≈e2,i for isotropic spatial distributions (corresponding to fruits), e0,i≥e1,i≈e2,i for predominantly linear distributions (branches), and e0,i≈e1,i≥e2,i for roughly planar distributions (leaves). In brief, the topological structure of PCD can be represented by ei=e0,ie1,ie2,i.


### 2.2.4. Tangent Space Character of PCD

Next, we design algorithm to calculate local tangent space character of PCD and determinek nearest neighbors pj′ of pi, pj′=p1′,p2′,…,pk′. We assume that a set of data points are sampled from a d-dimensional affine subspace; that is,(5)pj′=ci+Qiθj+εj1≤j≤k,where εj∈R3 represents noise vector, θj∈Rd is projection coordinates about pj′ on the local tangent space, and Θi=θ1,θ2,…,θk. d≤3. ci∈R3 is the origin coordinates of the tangent space and Qi∈R3×d is a matrix which forms an orthonormal basis of the affine subspace. The problem of linear manifold learning amounts to seek ci,Qi,θj to minimize the reconstruction error E; that is,(6)min⁡ci,Qi,θj∑j=1kpj-ci-Qiθj22=min⁡ci,Qi,ΘP-ci-QΘ22.The matrix of pi neighborhood is denoted by Xi=p1′,p2′,…,pk′, and we extract local information by calculating the eigenvectors and eigenvalue of the correlation matrix Xi-p-i′eTTXi-p-i′eT. Here p-i′=1/k∑j=1kpj, and e is a k-dimensional column vector of all ones; that is, (7)XiI-1k11TXiT=UiΛiUiT,where Ui  =ui1,ui2,…,uik is orthogonal matrix, and the diagonal elements of the diagonal matrix Λi are monotone decreasing, so the local tangent space information for the sample point pi is calculated:(8)ci=1kXi1,Qi=ui1,ui2,…,uid,Θi=QiTXiI-1k11T.From the above derivation, we can calculate the column vectoruim which is corresponding to the smallest diagonal element of Λi, and uim is also the normal vectors of local tangent space on pi.After the above analysis, a series of features about each pointpi are obtained, described as Fpi=v0,i,e0,i,e1,i,e2,i,uim, where v0,i=vi,x′,vi,y′,vi,z′ is the normal vector of pi, e0,i,e1,i,e2,i represent the distribution of v0,i, and uim are normal vectors of local tangent space on pi.


## 2.2.1. Data Collection by TLS

We used Leica Scan Station C10 to scan the same woods and obtain their PCD for processing. The specifications of C10 are described in Table1, and Figure 3 shows our experiment about tree data acquisition using Leica C10 scanner.Table 1
The specifications of Leica Scan Station C10.
C10 technical details
Accuracy of single measurement[position/distance/angular (horizontal/vertical)]
6 mm/4 mm/60μrad/60 μrad (12′′/12′′)
Range
300 m @ 90%; 134 m @ 18% albedo (minimum range 0.1 m)
Scan resolution
Spot size
From 0 to 50 m: 4.5 mm (FWHH-based);7 mm (Gaussian-based)
Point spacing
Fully selectable horizontal and vertical range; <1 mm minimum spacing
Laser plummet
Laser class: (IEC 60825-1)Centering accuracy: 1.5 mmLaser dotted diameter: 2.5 mm
Laser color
Green, wavelength = 532 nm visible
Environmental lighting
Fully operational between bright sunlight and complete darknessFigure 3
Using TLS to scan trees on our campus.
(a)
(b)The scanning PCD of woods are shown in Figure4. In order to calculate the effective leaf area index, we firstly design feature extraction methods to extract leaf PCD from the whole forestry PCD.Figure 4
The forestry PCD obtained by TLS.


## 2.2.2. PCD Features Extraction

Reliable feature extraction from 3D point cloud data is an important problem in many application domains, such as reverse engineering, object recognition, industrial inspection, and autonomous navigation. Many researchers have proposed many kinds of algorithms to extract the geometric features from 3D point cloud data [20, 25]. Tree’s PCD obtained by laser scanner have properties of both color and space location, but under the interference of external environment, such as illumination change and leaf swaying in the wind, the existing features of PCD, such as leaf, branch, and fruit, are inadequate to identify different tree organs. Therefore, in this paper new PCD features of topological structure and tangent space information are designed for PCD classification.


## 2.2.3. Topological Structure of PCD

We denote PCD set byP, P⊂R3, and pi=xi,yi,ziT is the point that belonged to set P. The k nearest points of pi are denoted by pj′=xj,yj,zjT with mean value u=1/k∑j=1kpj, and the covariance matrix of points pi and its neighbor is defined as Cp=1/N∑j=1Npj-upj-uT. After taking eigenvalue decomposition, such as Cpivi=λivi, the eigenvalue λi=λ0,i,λ1,i,λ2,i and eigenvector vi=v0,i,v1,i,v2,i of Cp are calculated, where the eigenvalue λi quantitatively shows the data variance along the axis vi or the compactness of the point distribution along the axis vi, so the v0,i corresponding to the smallest eigenvalue λ0,i is the normal vector of pi. Having got the norm vector of each PCD, we can use the following formula to calculate the covariance matrix Vpi of norm vector: (4)Vpi=1N∑j=1Nv0,j-v-iv0,j-v-iT,where v-i=1/N∑j=1Nvj is the average norm vector of pi’s neighbors. Consequently, we perform singular value decomposition on Vpi to identify the principal components of P’s spatial distribution. Specifically, we analyze the relative magnitudes of the top three eigenvalues ei=e0,ie1,ie2,i of Vpi. We expect e0,i≈e1,i≈e2,i for isotropic spatial distributions (corresponding to fruits), e0,i≥e1,i≈e2,i for predominantly linear distributions (branches), and e0,i≈e1,i≥e2,i for roughly planar distributions (leaves). In brief, the topological structure of PCD can be represented by ei=e0,ie1,ie2,i.


## 2.2.4. Tangent Space Character of PCD

Next, we design algorithm to calculate local tangent space character of PCD and determinek nearest neighbors pj′ of pi, pj′=p1′,p2′,…,pk′. We assume that a set of data points are sampled from a d-dimensional affine subspace; that is,(5)pj′=ci+Qiθj+εj1≤j≤k,where εj∈R3 represents noise vector, θj∈Rd is projection coordinates about pj′ on the local tangent space, and Θi=θ1,θ2,…,θk. d≤3. ci∈R3 is the origin coordinates of the tangent space and Qi∈R3×d is a matrix which forms an orthonormal basis of the affine subspace. The problem of linear manifold learning amounts to seek ci,Qi,θj to minimize the reconstruction error E; that is,(6)min⁡ci,Qi,θj∑j=1kpj-ci-Qiθj22=min⁡ci,Qi,ΘP-ci-QΘ22.The matrix of pi neighborhood is denoted by Xi=p1′,p2′,…,pk′, and we extract local information by calculating the eigenvectors and eigenvalue of the correlation matrix Xi-p-i′eTTXi-p-i′eT. Here p-i′=1/k∑j=1kpj, and e is a k-dimensional column vector of all ones; that is, (7)XiI-1k11TXiT=UiΛiUiT,where Ui  =ui1,ui2,…,uik is orthogonal matrix, and the diagonal elements of the diagonal matrix Λi are monotone decreasing, so the local tangent space information for the sample point pi is calculated:(8)ci=1kXi1,Qi=ui1,ui2,…,uid,Θi=QiTXiI-1k11T.From the above derivation, we can calculate the column vectoruim which is corresponding to the smallest diagonal element of Λi, and uim is also the normal vectors of local tangent space on pi.After the above analysis, a series of features about each pointpi are obtained, described as Fpi=v0,i,e0,i,e1,i,e2,i,uim, where v0,i=vi,x′,vi,y′,vi,z′ is the normal vector of pi, e0,i,e1,i,e2,i represent the distribution of v0,i, and uim are normal vectors of local tangent space on pi.


## 2.3. PCD Classification Using GMM and EM Algorithm

The scanning PCD data belong to the two semantic classes (branch and leaf). In this section, we combine Gaussian mixture model (GMM) classifier and expectation maximization (EM) algorithm to classify forestry PCD.GMM is a weighted sum ofA component Gaussian densities as given by (9),(9)pFp∣λ=∑i=1AωigFp∣μi,σi,where Fp is a continuously valued data vector of PCD (i.e., measurement or features), ωi,i=1,…,A, are the mixture weights, and gFp∣μi,σi, i=1,…,A, are the component Gaussian densities. PCD can be divided into two categories, which are branch and leaf, respectively, so A=2. Each component density is a form of Gaussian function, with mean vector μi and covariance matrix σi, and ωi is the weight coefficient of each class. The expansion formula of g is defined as follows: (10)gFp∣μi,σi=12πσi1/2·exp⁡-12Fp-μi′σi-1Fp-μi.The mixture weights satisfy the constraint that∑i=1Mωi=1. The complete Gaussian mixture model is parameterized by the mean vectors, covariance matrices, and mixture weights from all component densities. These parameters are collectively represented by the notation, λ=ωi,μi,σi,i=1,…,A. Then the expectation maximization (EM) algorithm is proposed to maximize the likelihood pFp∣λ of the data Fp drawn from an unknown distribution. In each iterative process of the algorithm, two steps of EM are executed. Firstly, estimate the distribution of the hidden variable according to the current value of the parameters; secondly, modify the parameters in order to maximize the joint distribution of the data and the hidden variable. Specific formula is expressed as follows:(11)λ∗=argmax⁡λ∏j=1N′∑k=1AωkgFp∣λ,where N′ represents total number of the calculated PCD. The output of the GMM and EM generates an independent classification for each of the 3D PCD based solely upon PCD’s feature vectors and confirmed in our experiments; such classification can be quite effective and the different plant organs are distinguished; thus, the whole PCD can be divided into two types: branches Pbranch and leaves Pleaf.


## 2.4. Spherical Projection of Classified PCD

In this section, we convert the segmented woods PCD from Cartesian coordinate system to spherical coordinate system with a radius of one, which projects the PCD onto the surface of a hemisphere. The aim of the stereographic projection was to project all points from the surface of the upper hemisphere (3D space) into the projection plane that is a 2D space. For example, the PCDPleafxi,yi,zi projected onto the spherical surface would be Pleaf′ri,φi,θi, where φ is azimuth, and θ is zenith angle. The projection process is computed and converted using (12)ri=xi2+yi2+zi2,φi=arctanyixi,θi=arccosziri.After the projection transform, each pointPleaf′ri,φi,θi is projected onto the X-Y plane and the hemispherical photograph Ifisheye′ transformed from PCD is obtained. In the following section, we use forestry models to calculate canopy indicators.


## 3. Canopy Indicator Retrieval

### 3.1. Gap Fraction and Clumping Index Calculation

According to formula (1) in Section 2.1, the gap fraction Tθ of Ifisheye and Ifisheye′ under different viewing zenith angle is calculated:(13)Tθ=sumBi,j,θ==0sumBi,j,θ~=111,B∈IfisheyeorIfisheye′.Many optical instruments measure canopy gap fraction based on radiation transmission through the canopy. Assuming random spatial distribution of leaf, the effective LA1 can be calculated from the gap fraction by adopting Miller’s theorem [26] summarized in the following equation:(14)LAI=2∫0π/2ln⁡1Tθcos⁡θsinθdθ.Discrete integral form of formula (14) is shown as follows:(15)LAI=2∑i=1nln⁡1Tθicos⁡θisinθiΔθ,where n represented the number of rings which divided the HP into small sections according to the zenith angle and is mentioned in Section 2.1. Consequently, the clumping index Ω can be defined by logarithmic mean value of the gap fraction, and specific formula is expressed as follows:(16)Ωθ=ln⁡Tθ,φ¯ln⁡Tθ,φ¯.


### 3.2. Leaf Inclination Angle Calculation Based on PCD

The branch and leaf PCD, denoted byPbranch and Pleaf, are classified by the method of Section 2. In addition, the normal vector v0,i=vi,x′,vi,y′,vi,z′ of each leaf point pi∈Pleaf is obtained, so the leaf inclination angle αi can be described as the included angle between v0,i and Z⃑=0,0,1 axis, which is defined by formula (17) and shown in Figure 5:(17)αi=arccosvi,x′·0+vi,y′·0+vi,z′·1vi,x′2+vi,y′2+vi,z′2∗02+02+12.The average leaf inclination angle α- can be found by computing the following equation:(18)α-=∑i=1PleafαiPleaf,where Pleaf is the number of points in the cloud.Figure 5
Leaf inclination angle calculation based on PCD.


### 3.3. LAI Calculation Based on Ellipsoid Model

Leaf oriented distribution is a critical component for estimating LAI from optical measurements. The LAI estimation process is usually computed based on the following equation which derived from Beer-Lambert law:(19)LAI=-ln⁡Tθcos⁡θGθ,α-,where α- represent average leaf inclination angle, where Tθ is the gap fraction defined as the probability of abeam transmitted through the canopy with an incident angle θ.The following ellipsoidal distribution function [27] is chosen for Gθ,α- in this study for its simplicity and reasonable generality. Gθ,α- is the extinction coefficient, defined as the mean projection of unit foliage area on the plane normal to the direction of the beam and specified as formula (20):(20)Gθ,α-=x2+tan2θ1/2cos⁡θx+sin-1ε1/ε1,x≤1,Gθ,α-=x2+tan2θ1/2cos⁡θx+1/2ε2xln⁡1+ε2/1-ε2,x>1,where ε1=1-x21/2 and ε2=1-x-21/2. The value of x is related to the vertical a and horizontal b axes of the ellipsoid by x=b/a. When x=1 the distribution is spherical (random) and Gθ,α-=sinθ. As x increases or decreases from 1, the area distribution becomes planophile or erectophile, respectively [28]. In addition, we used ellipsoidal leaf angle distribution [29] to describe the relationship between the leaf orientation and the zenith angle of the radiation. x can be estimated by average leaf orientation: (21)x=-3+α-9.65-0.6061.If a set ofα- values are assumed to be α-=5∘,15∘,25∘,35∘,…, we can draw the diagram illustrating the Gθ,α- values versus the variation of θ as shown in Figure 6. Observed Gθ,α- appear to be a relatively stable value of 0.5 when θ≈57∘, so the LAI can be calculated by the following formula as θ is equal to 57.5∘:(22)LAI=-ln⁡T57.5∘cos⁡57.5∘G57.5∘.Figure 6
Extinction coefficient versus average leaf inclination angle and beam incident angle.As mentioned earlier, the average leaf inclination angleα- can be obtained by (18). Meanwhile, x value can be calculated by using (21), and it follows that the ellipsoidal distribution of leaf angles can be represented as follows: (23)Pα=2x3sinασcos⁡2α+x2sin2α2,where α is the leaf inclination angle and (24)σ=x+arcsin1-x21-x2,x<1,σ=x+ln⁡1+1-x-2/1-1-x-22x1-x-2,x≥1,x=-3+α-9.65-0.6061.Using (23), we can calculate PCD angles distribution of different tree organs, which can be compared with our method of PCD features extraction. The similar results have been obtained as shown in Figures 8(c) and 8(h).


### 3.4. LAI Retrieval by Iterative Inversion Method

By substituting (20) and (21) into (19), we can get (25)Tsimθ=eτ,τ=LAIx2+tan2θx+(sin-11-x2/1-x2),x=-3+α-9.65-0.6061.Supposing the initial values of LAI,α- and θ are known, and the simulative gap fraction Tsim(θ) can be obtained by (25). Meanwhile, we have got canopy gap fraction Tθ of each solar zenith angle θ using formula (3). Consequently, we adjust the value of LAI and α- to satisfy the minimization of (26) and thus the appropriate quantity of LAI and α- that can be obtained as the ideal results of the forestry indicators:(26)LA^I,α-^=argmin⁡A,α-∑θTθ-Tsimθ2.


### 3.5. LAI Estimation Based on the Clumping Index

Another expression form of formula (19) can be written as (27)LAIt·Ω=LAI=-ln⁡Tθcos⁡θGθ,α-,where LAIt is the true leaf area index and Ω is the clumping index [30]. Specific formula is as follows:(28)LAIt=LAIΩ,(29)Ωθ=ln⁡Tθ,φ¯ln⁡Tθ,φ¯.Through calculating the logarithmic mean of gap fraction, we can get clumping indexΩ and true leaf area index LAIt.


## 3.1. Gap Fraction and Clumping Index Calculation

According to formula (1) in Section 2.1, the gap fraction Tθ of Ifisheye and Ifisheye′ under different viewing zenith angle is calculated:(13)Tθ=sumBi,j,θ==0sumBi,j,θ~=111,B∈IfisheyeorIfisheye′.Many optical instruments measure canopy gap fraction based on radiation transmission through the canopy. Assuming random spatial distribution of leaf, the effective LA1 can be calculated from the gap fraction by adopting Miller’s theorem [26] summarized in the following equation:(14)LAI=2∫0π/2ln⁡1Tθcos⁡θsinθdθ.Discrete integral form of formula (14) is shown as follows:(15)LAI=2∑i=1nln⁡1Tθicos⁡θisinθiΔθ,where n represented the number of rings which divided the HP into small sections according to the zenith angle and is mentioned in Section 2.1. Consequently, the clumping index Ω can be defined by logarithmic mean value of the gap fraction, and specific formula is expressed as follows:(16)Ωθ=ln⁡Tθ,φ¯ln⁡Tθ,φ¯.


## 3.2. Leaf Inclination Angle Calculation Based on PCD

The branch and leaf PCD, denoted byPbranch and Pleaf, are classified by the method of Section 2. In addition, the normal vector v0,i=vi,x′,vi,y′,vi,z′ of each leaf point pi∈Pleaf is obtained, so the leaf inclination angle αi can be described as the included angle between v0,i and Z⃑=0,0,1 axis, which is defined by formula (17) and shown in Figure 5:(17)αi=arccosvi,x′·0+vi,y′·0+vi,z′·1vi,x′2+vi,y′2+vi,z′2∗02+02+12.The average leaf inclination angle α- can be found by computing the following equation:(18)α-=∑i=1PleafαiPleaf,where Pleaf is the number of points in the cloud.Figure 5
Leaf inclination angle calculation based on PCD.


## 3.3. LAI Calculation Based on Ellipsoid Model

Leaf oriented distribution is a critical component for estimating LAI from optical measurements. The LAI estimation process is usually computed based on the following equation which derived from Beer-Lambert law:(19)LAI=-ln⁡Tθcos⁡θGθ,α-,where α- represent average leaf inclination angle, where Tθ is the gap fraction defined as the probability of abeam transmitted through the canopy with an incident angle θ.The following ellipsoidal distribution function [27] is chosen for Gθ,α- in this study for its simplicity and reasonable generality. Gθ,α- is the extinction coefficient, defined as the mean projection of unit foliage area on the plane normal to the direction of the beam and specified as formula (20):(20)Gθ,α-=x2+tan2θ1/2cos⁡θx+sin-1ε1/ε1,x≤1,Gθ,α-=x2+tan2θ1/2cos⁡θx+1/2ε2xln⁡1+ε2/1-ε2,x>1,where ε1=1-x21/2 and ε2=1-x-21/2. The value of x is related to the vertical a and horizontal b axes of the ellipsoid by x=b/a. When x=1 the distribution is spherical (random) and Gθ,α-=sinθ. As x increases or decreases from 1, the area distribution becomes planophile or erectophile, respectively [28]. In addition, we used ellipsoidal leaf angle distribution [29] to describe the relationship between the leaf orientation and the zenith angle of the radiation. x can be estimated by average leaf orientation: (21)x=-3+α-9.65-0.6061.If a set ofα- values are assumed to be α-=5∘,15∘,25∘,35∘,…, we can draw the diagram illustrating the Gθ,α- values versus the variation of θ as shown in Figure 6. Observed Gθ,α- appear to be a relatively stable value of 0.5 when θ≈57∘, so the LAI can be calculated by the following formula as θ is equal to 57.5∘:(22)LAI=-ln⁡T57.5∘cos⁡57.5∘G57.5∘.Figure 6
Extinction coefficient versus average leaf inclination angle and beam incident angle.As mentioned earlier, the average leaf inclination angleα- can be obtained by (18). Meanwhile, x value can be calculated by using (21), and it follows that the ellipsoidal distribution of leaf angles can be represented as follows: (23)Pα=2x3sinασcos⁡2α+x2sin2α2,where α is the leaf inclination angle and (24)σ=x+arcsin1-x21-x2,x<1,σ=x+ln⁡1+1-x-2/1-1-x-22x1-x-2,x≥1,x=-3+α-9.65-0.6061.Using (23), we can calculate PCD angles distribution of different tree organs, which can be compared with our method of PCD features extraction. The similar results have been obtained as shown in Figures 8(c) and 8(h).


## 3.4. LAI Retrieval by Iterative Inversion Method

By substituting (20) and (21) into (19), we can get (25)Tsimθ=eτ,τ=LAIx2+tan2θx+(sin-11-x2/1-x2),x=-3+α-9.65-0.6061.Supposing the initial values of LAI,α- and θ are known, and the simulative gap fraction Tsim(θ) can be obtained by (25). Meanwhile, we have got canopy gap fraction Tθ of each solar zenith angle θ using formula (3). Consequently, we adjust the value of LAI and α- to satisfy the minimization of (26) and thus the appropriate quantity of LAI and α- that can be obtained as the ideal results of the forestry indicators:(26)LA^I,α-^=argmin⁡A,α-∑θTθ-Tsimθ2.


## 3.5. LAI Estimation Based on the Clumping Index

Another expression form of formula (19) can be written as (27)LAIt·Ω=LAI=-ln⁡Tθcos⁡θGθ,α-,where LAIt is the true leaf area index and Ω is the clumping index [30]. Specific formula is as follows:(28)LAIt=LAIΩ,(29)Ωθ=ln⁡Tθ,φ¯ln⁡Tθ,φ¯.Through calculating the logarithmic mean of gap fraction, we can get clumping indexΩ and true leaf area index LAIt.


## 4. Experimental Results

### 4.1. Experimental Setup

There is a dense forest behind Nanjing Forestry University, and we randomly chose stand sample plots in this forestry as the experimental subject. Then, the experiment was set up at this sample plot with some tree species such asPterocarya stenoptera C. DC, Sect.Leucoides Spach, andSophora japonica Linn. We, respectively, used three devices such as TLS, fisheye camera, and LAI-2200 plant canopy analyzer [31] to scan the same trees at fixed position and fixed angle of 90 degrees, just as shown in Figure 7(a). Forestry field measurements were performed sequentially for the many plots to capture the canopy structural dynamics. All instruments were newly procured and factory calibrated. All optical measurements were conducted near sunset or under overcast conditions because the parameter sensitivity and the retrieval errors increase under direct illuminations. Meanwhile, we adopted our methods to compute leaf area index, mean tilt angle, canopy gap fraction, and so forth from obtained data and compared all the results.Figure 7
Calculation of canopy gap fraction from hemispherical photographs.
(a)
(b)
(c)
(d)
(e)Figure 8
Canopy parameter calculation based on the laser PCD.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)


### 4.2. Experiment Using Fisher Eye Camera

We took hemispherical photographs of the forestry canopy by fisheye lens and Canon 5d (Mark III) camera and used our algorithm to process the image in order to calculate the relevant forestry indexes. The specific details are as follows.In Figure7, (b) is the hemispherical photograph taken by our devices. After our image processing algorithm, the result is shown in (c) that branch and background pixels are removed. The gap fraction versus different zenith angle, illustrated by red lines in (d), is calculated from 12 different hemispherical photographs of the same forestry canopy. Meanwhile, we also used LAI-2200 to measure the gap fraction at 5 zenith angles such as 7∘, 23∘, 38∘, 53∘, and 68∘, which basically coincided with the results of HP and is shown by blue lines. If we change the number of concentric rings of the HP, the value of variable n changes from 9 to 288, and the corresponding gap fraction results of the same HP are shown in (e).


### 4.3. Experiment Using TLS

This trial used TLS to capture PCD of forest stand and adopted classification and projection algorithm as originally proposed by us to calculate canopy indexes. In addition, we adopted LAI-2200 to calculate parameters of the same forest stand and compared the results of various sampling plots obtained by these devices. The specific details are illustrated in Figure8.Figure8(a) shows the scanned PCD of the same forest plot. (b) shows the classification result of our algorithm that divides the PCD into two class: leaf and branch. (c) is leaf angle distribution diagram, in which yellow and green histograms, respectively, represent leaf and branch angle distribution of dense PCD, and blue and red histograms, respectively, represent leaf and branch angle distribution of sparse PCD. The inclination angle of leaf PCD presents uniform distribution, but branch point shows vertical properties, so the number of branch points in cloud gradually increased as the angle approaches 90 degrees. The forestry PCD after classification is projected onto the hemisphere coordinates (Figure 8(d)) and hemisphere plane (Figure 8(e)). By using image segmentation method to process Figure 8(e), we can get the corresponding gap fraction values versus the variable of zenith angle, which is shown in (f). Comparison diagram (g) depicts calculated gap fraction versus different zenith angel with three methods: hemispherical photography, terrestrial laser scanning, and LAI-2200 plant canopy analyzer. While zenith angle varies at the range from 30 to 90 degrees (labeled by green dotted lines and most relevant to LAI estimation), we can find that the gap fraction obtained by these three methods got the same results. The calculation results of formula (23) are displayed by (h), which depicts inclination angle distribution of branch and leaf PCD, and the shapes of these curves are similar to the histogram distribution in (c). By substituting (3) and (20) into (19), (i) illustrates the result of LAI values obtained from PCD. In all six PCD types, a range of high, medium, and low densities was recorded, which reflects the variation of LAI estimation using our methods. When the zenith angle varies within the range from 40 to 70 degrees, the LAI value estimated by terrestrial laser scanning method tended to be similar to the result of hemispherical photography method and LAI-2200 device, as labeled by red dotted lines and red pentagram in (i), respectively. In order to reflect the statistical testing of our work, we adopted sampling strategy that allows us to gain information about a site without looking at every plot in the forestry. Sampling involves taking measurements on some small plots that are representative of the larger study area and using the data to represent the entire study area. By locating quadrats using random or stratified sampling, we simultaneously chose other plots in the forestry and used the three methods to obtain each quadrat’s index. Comparing with LAI retrieved from different optical methods, we drew Figure 8(j) to show a very good relationship among the LAI values estimated using HP and PCD methods (R2=0.7052, RMSE = 0.1323). A good relationship between the LAI-2200 5-ring and PCD (R2=0.5742, RMSE = 0.1590) can also be observed in Figure 8(k).


### 4.4. Experimental Result Analysis

Using our three methods, hemispherical photography, LAI-2200 plant canopy analyzer, and terrestrial laser scanning, the experimental results are obtained and summarized in Tables2–5. In Tables 2 and 5, the total parameters of the same forest stand are estimated from HP and PCD projection image, combining with theoretical approaches mentioned above, such as Miller’s method (equation (14)), Iterative inversion method (equation (26)), LAI calculation according to single zenith angle (equation (22)), and canopy clumping index calculation by (29). In addition, we examine the effect of our results on the calculation accuracy with LAI-2200 plant canopy analyzer (Table 4). Judging by these tables and Figures 8(j) and 8(k), it can be found that the three methods get similar results of measuring forestry stand, as labeled by digits in bold. Conversely, we also find that each method has its own disadvantages. Firstly, the hemispherical photography based gap fraction estimation method will be affected by the light environment in the forest stands and camera exposure index. For example, cloudy skies are ideal light conditions for digital hemispherical photograph to differentiate foliage elements from the sky background and branch elements of hemispherical photographs. Considering the effectiveness of LAI-2200 device, one of the traditional underlying assumptions has been that foliage absorbs all the radiation in the blue waveband seen by the sensor (320–490 nm). This is usually a good assumption under diffuse light conditions such as uniform overcast, just before sunrise, or just after sunset. In direct sunlight, however, reflectance of foliage causes a much greater overestimation of the gap fraction and underestimation of leaf area index, so the mechanism for correcting measurements for the radiation reflected and transmitted by the foliage is necessarily dependent when the scattering errors are the highest. In short, results measured by LAI-2200 still have deviation and depend on artificial correction, whereas, the TLS-based gap fraction estimation method using laser can work in both cloudy and sunny light conditions. The measurement errors of TLS are affected by the density of the sampling points and size of laser beam spot, as well as the accuracy of classification algorithm for leaf and branch PCD recognition. Secondly, combining with salient feature extraction from PCD, the PCD can be applied not only to estimate the gap fraction by spherical projection method but also to estimate leaf inclination angle from normal vector calculation of PCD. In addition, terrestrial laser scanning method can estimate the LAI from a 3D perspective and permanently record the 3D structural information of forest canopy. Thirdly, the abundant information implicitly contained within the PCD generated using TLS could provide more information about the canopy structure besides LAI. This includes timber volume, leaf area density, leaf morphology, carbon content, and canopy volume. The additional information will enhance our ability to monitor forest tree structure from a 3D perspective dynamically, which could provide useful data for the long-term ecological studies.Table 2
LAI calculation based on hemispherical photography.
Data
Gap fraction of 57° degrees
LAI
Canopy clumping index
Miller’s method
Iterative inversion method
Equation (22) of singe angle
HP 1
0.2132
1.4893
1.7204
1.6606
0.9338
HP 2
0.1802
1.4924
1.7261
1.8416
0.9375
HP 3
0.1746
1.4813
1.7067
1.8755
0.9406
HP 4
0.1990
1.4860
1.6993
1.7346
0.9304
HP 5
0.1954
1.4828
1.7158
1.7544
0.9277
HP 6
0.1970
1.4822
1.7099
1.7460
0.9269
HP 7
0.1696
1.4493
1.6663
1.9070
0.9412
HP 8
0.1647
1.4716
1.6976
1.9382
0.9354
HP 9
0.1699
1.4623
1.6744
1.9047
0.9374
HP 10
0.1694
1.4829
1.6998
1.9079
0.9319
HP 11
0.1741
1.4876
1.7117
1.8786
0.9295
HP 12
0.1699
1.4696
1.6992
1.9046
0.9343Table 3
Average leaf inclination angle calculation based on PCD.
Experimental data
Average inclination angle of branch PCD
Average inclination angle of leaf PCD
Medium-density PCD1
59.5360 degrees
54.5971 degrees
High-density PCD1
58.2641 degrees
54.2190 degrees
Low-density PCD1
67.9814 degrees
53.4913 degrees
Medium-density PCD2
58.4761 degrees
54.1445 degrees
Low-density PCD2
58.7912 degrees
54.1732 degrees
High-density PCD2
69.5972 degrees
53.8981 degrees
Average inclination angle
62.1077 degrees
54.0872 degreesTable 4
Parameters obtained by LAI-2200.
Case
Mean leaf inclination angle
Diffusion noninterception
LAI and standard deviation
1
49.536 degrees
0.184
1.938 ± 0.22
2
52.732 degrees
0.233
2.105 ± 0.14
3
51.130 degrees
0.211
1.921 ± 0.19
4
52.178 degrees
0.194
2.224 ± 0.20
5
50.791 degrees
0.207
1.883 ± 0.17
6
49.870 degrees
0.193
2.001 ± 0.12Table 5
LAI calculation based on PCD.
PCD data
Gap fraction of 57 degrees
LAI
Canopy clumping index
Miller’s method
Iterative  inversion method
Equation  (22) of singe angle
Low-density PCD1
0.3460
1.2169
1.2887
1.1661
0.9778
Medium-density PCD1
0.1393
2.2692
2.4095
2.2092
0.9296
High-density PCD1
0.1283
2.4099
2.6309
2.3559
0.9138
Low-density PCD2
0.3438
1.2224
1.3383
1.1662
0.9761
Medium-density PCD2
0.1894
1.9356
2.1179
1.8581
0.9557
High-density PCD2
0.1184
2.4654
2.6572
2.3485
0.9328


## 4.1. Experimental Setup

There is a dense forest behind Nanjing Forestry University, and we randomly chose stand sample plots in this forestry as the experimental subject. Then, the experiment was set up at this sample plot with some tree species such asPterocarya stenoptera C. DC, Sect.Leucoides Spach, andSophora japonica Linn. We, respectively, used three devices such as TLS, fisheye camera, and LAI-2200 plant canopy analyzer [31] to scan the same trees at fixed position and fixed angle of 90 degrees, just as shown in Figure 7(a). Forestry field measurements were performed sequentially for the many plots to capture the canopy structural dynamics. All instruments were newly procured and factory calibrated. All optical measurements were conducted near sunset or under overcast conditions because the parameter sensitivity and the retrieval errors increase under direct illuminations. Meanwhile, we adopted our methods to compute leaf area index, mean tilt angle, canopy gap fraction, and so forth from obtained data and compared all the results.Figure 7
Calculation of canopy gap fraction from hemispherical photographs.
(a)
(b)
(c)
(d)
(e)Figure 8
Canopy parameter calculation based on the laser PCD.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)


## 4.2. Experiment Using Fisher Eye Camera

We took hemispherical photographs of the forestry canopy by fisheye lens and Canon 5d (Mark III) camera and used our algorithm to process the image in order to calculate the relevant forestry indexes. The specific details are as follows.In Figure7, (b) is the hemispherical photograph taken by our devices. After our image processing algorithm, the result is shown in (c) that branch and background pixels are removed. The gap fraction versus different zenith angle, illustrated by red lines in (d), is calculated from 12 different hemispherical photographs of the same forestry canopy. Meanwhile, we also used LAI-2200 to measure the gap fraction at 5 zenith angles such as 7∘, 23∘, 38∘, 53∘, and 68∘, which basically coincided with the results of HP and is shown by blue lines. If we change the number of concentric rings of the HP, the value of variable n changes from 9 to 288, and the corresponding gap fraction results of the same HP are shown in (e).


## 4.3. Experiment Using TLS

This trial used TLS to capture PCD of forest stand and adopted classification and projection algorithm as originally proposed by us to calculate canopy indexes. In addition, we adopted LAI-2200 to calculate parameters of the same forest stand and compared the results of various sampling plots obtained by these devices. The specific details are illustrated in Figure8.Figure8(a) shows the scanned PCD of the same forest plot. (b) shows the classification result of our algorithm that divides the PCD into two class: leaf and branch. (c) is leaf angle distribution diagram, in which yellow and green histograms, respectively, represent leaf and branch angle distribution of dense PCD, and blue and red histograms, respectively, represent leaf and branch angle distribution of sparse PCD. The inclination angle of leaf PCD presents uniform distribution, but branch point shows vertical properties, so the number of branch points in cloud gradually increased as the angle approaches 90 degrees. The forestry PCD after classification is projected onto the hemisphere coordinates (Figure 8(d)) and hemisphere plane (Figure 8(e)). By using image segmentation method to process Figure 8(e), we can get the corresponding gap fraction values versus the variable of zenith angle, which is shown in (f). Comparison diagram (g) depicts calculated gap fraction versus different zenith angel with three methods: hemispherical photography, terrestrial laser scanning, and LAI-2200 plant canopy analyzer. While zenith angle varies at the range from 30 to 90 degrees (labeled by green dotted lines and most relevant to LAI estimation), we can find that the gap fraction obtained by these three methods got the same results. The calculation results of formula (23) are displayed by (h), which depicts inclination angle distribution of branch and leaf PCD, and the shapes of these curves are similar to the histogram distribution in (c). By substituting (3) and (20) into (19), (i) illustrates the result of LAI values obtained from PCD. In all six PCD types, a range of high, medium, and low densities was recorded, which reflects the variation of LAI estimation using our methods. When the zenith angle varies within the range from 40 to 70 degrees, the LAI value estimated by terrestrial laser scanning method tended to be similar to the result of hemispherical photography method and LAI-2200 device, as labeled by red dotted lines and red pentagram in (i), respectively. In order to reflect the statistical testing of our work, we adopted sampling strategy that allows us to gain information about a site without looking at every plot in the forestry. Sampling involves taking measurements on some small plots that are representative of the larger study area and using the data to represent the entire study area. By locating quadrats using random or stratified sampling, we simultaneously chose other plots in the forestry and used the three methods to obtain each quadrat’s index. Comparing with LAI retrieved from different optical methods, we drew Figure 8(j) to show a very good relationship among the LAI values estimated using HP and PCD methods (R2=0.7052, RMSE = 0.1323). A good relationship between the LAI-2200 5-ring and PCD (R2=0.5742, RMSE = 0.1590) can also be observed in Figure 8(k).


## 4.4. Experimental Result Analysis

Using our three methods, hemispherical photography, LAI-2200 plant canopy analyzer, and terrestrial laser scanning, the experimental results are obtained and summarized in Tables2–5. In Tables 2 and 5, the total parameters of the same forest stand are estimated from HP and PCD projection image, combining with theoretical approaches mentioned above, such as Miller’s method (equation (14)), Iterative inversion method (equation (26)), LAI calculation according to single zenith angle (equation (22)), and canopy clumping index calculation by (29). In addition, we examine the effect of our results on the calculation accuracy with LAI-2200 plant canopy analyzer (Table 4). Judging by these tables and Figures 8(j) and 8(k), it can be found that the three methods get similar results of measuring forestry stand, as labeled by digits in bold. Conversely, we also find that each method has its own disadvantages. Firstly, the hemispherical photography based gap fraction estimation method will be affected by the light environment in the forest stands and camera exposure index. For example, cloudy skies are ideal light conditions for digital hemispherical photograph to differentiate foliage elements from the sky background and branch elements of hemispherical photographs. Considering the effectiveness of LAI-2200 device, one of the traditional underlying assumptions has been that foliage absorbs all the radiation in the blue waveband seen by the sensor (320–490 nm). This is usually a good assumption under diffuse light conditions such as uniform overcast, just before sunrise, or just after sunset. In direct sunlight, however, reflectance of foliage causes a much greater overestimation of the gap fraction and underestimation of leaf area index, so the mechanism for correcting measurements for the radiation reflected and transmitted by the foliage is necessarily dependent when the scattering errors are the highest. In short, results measured by LAI-2200 still have deviation and depend on artificial correction, whereas, the TLS-based gap fraction estimation method using laser can work in both cloudy and sunny light conditions. The measurement errors of TLS are affected by the density of the sampling points and size of laser beam spot, as well as the accuracy of classification algorithm for leaf and branch PCD recognition. Secondly, combining with salient feature extraction from PCD, the PCD can be applied not only to estimate the gap fraction by spherical projection method but also to estimate leaf inclination angle from normal vector calculation of PCD. In addition, terrestrial laser scanning method can estimate the LAI from a 3D perspective and permanently record the 3D structural information of forest canopy. Thirdly, the abundant information implicitly contained within the PCD generated using TLS could provide more information about the canopy structure besides LAI. This includes timber volume, leaf area density, leaf morphology, carbon content, and canopy volume. The additional information will enhance our ability to monitor forest tree structure from a 3D perspective dynamically, which could provide useful data for the long-term ecological studies.Table 2
LAI calculation based on hemispherical photography.
Data
Gap fraction of 57° degrees
LAI
Canopy clumping index
Miller’s method
Iterative inversion method
Equation (22) of singe angle
HP 1
0.2132
1.4893
1.7204
1.6606
0.9338
HP 2
0.1802
1.4924
1.7261
1.8416
0.9375
HP 3
0.1746
1.4813
1.7067
1.8755
0.9406
HP 4
0.1990
1.4860
1.6993
1.7346
0.9304
HP 5
0.1954
1.4828
1.7158
1.7544
0.9277
HP 6
0.1970
1.4822
1.7099
1.7460
0.9269
HP 7
0.1696
1.4493
1.6663
1.9070
0.9412
HP 8
0.1647
1.4716
1.6976
1.9382
0.9354
HP 9
0.1699
1.4623
1.6744
1.9047
0.9374
HP 10
0.1694
1.4829
1.6998
1.9079
0.9319
HP 11
0.1741
1.4876
1.7117
1.8786
0.9295
HP 12
0.1699
1.4696
1.6992
1.9046
0.9343Table 3
Average leaf inclination angle calculation based on PCD.
Experimental data
Average inclination angle of branch PCD
Average inclination angle of leaf PCD
Medium-density PCD1
59.5360 degrees
54.5971 degrees
High-density PCD1
58.2641 degrees
54.2190 degrees
Low-density PCD1
67.9814 degrees
53.4913 degrees
Medium-density PCD2
58.4761 degrees
54.1445 degrees
Low-density PCD2
58.7912 degrees
54.1732 degrees
High-density PCD2
69.5972 degrees
53.8981 degrees
Average inclination angle
62.1077 degrees
54.0872 degreesTable 4
Parameters obtained by LAI-2200.
Case
Mean leaf inclination angle
Diffusion noninterception
LAI and standard deviation
1
49.536 degrees
0.184
1.938 ± 0.22
2
52.732 degrees
0.233
2.105 ± 0.14
3
51.130 degrees
0.211
1.921 ± 0.19
4
52.178 degrees
0.194
2.224 ± 0.20
5
50.791 degrees
0.207
1.883 ± 0.17
6
49.870 degrees
0.193
2.001 ± 0.12Table 5
LAI calculation based on PCD.
PCD data
Gap fraction of 57 degrees
LAI
Canopy clumping index
Miller’s method
Iterative  inversion method
Equation  (22) of singe angle
Low-density PCD1
0.3460
1.2169
1.2887
1.1661
0.9778
Medium-density PCD1
0.1393
2.2692
2.4095
2.2092
0.9296
High-density PCD1
0.1283
2.4099
2.6309
2.3559
0.9138
Low-density PCD2
0.3438
1.2224
1.3383
1.1662
0.9761
Medium-density PCD2
0.1894
1.9356
2.1179
1.8581
0.9557
High-density PCD2
0.1184
2.4654
2.6572
2.3485
0.9328


## 5. Conclusions and Future Work

In order to nondestructively estimate the forestry indicators of heterogeneous canopies, we used three methods, hemispherical photography, LAI-2200 plant canopy analyzer, and terrestrial laser scanning, to acquire high dimensional data of forestry stand from different perspectives. In addition, combining with forestry empirical models, we design geometrical projection methods, PCD classification, and image segmentation algorithm to estimate parameters of forestry stand. Our results have shown that the TLS could be used to estimate LAI for heterogeneous forests at forest plot level without the restriction of the light environment. Meanwhile, we are able to demonstrate relationships between HP, LAI-2200, and TLS for canopy index estimation. The relationship also suggests that TLS could serve as a calibration tool for HP-based and airborne LiDAR-based LAI estimation. In the future work, we will adopt LiDAR, high spectral radiometer and tracing radiation canopy architecture analyzer, and TLS to comprehensively obtain forestry stand data with forestry-related calculation model proposed to estimate forestry parameters. As further development of these tools and methods for extracting biophysical and ecological parameters from TLS data sets, long-term forest ecosystem monitoring will benefit from repeatable techniques assuring data for sustainable forest management practices.



---

*Source: 206108-2015-10-05.xml*