Comparative Evaluation of Algorithms for Leaf Area Index Estimation from Digital Hemispherical Photography through Virtual Forests

The in situ leaf area index (LAI) measurement plays a vital role in calibrating and validating satellite LAI products. Digital hemispherical photography (DHP) is a widely used in situ forest LAI measurement method. There have been many software programs encompassing a variety of algorithms to estimate LAI from DHP. However, there is no conclusive study for an accuracy comparison among them, due to the difficulty in acquiring forest LAI reference values. In this study, we aim to use virtual (i.e., computer-simulated) broadleaf forests for the accuracy assessment of LAI algorithms in commonly used LAI software programs. Three commonly used DHP programs, including Can_Eye, CIMES, and Hemisfer, were selected since they provide estimates of both effective LAI and true LAI. Individual tree models with and without leaves were first reconstructed based on terrestrial LiDAR point clouds. Various stands were then created from these models. A ray-tracing technique was combined with the virtual forests to model synthetic DHP, for both leaf-on and leaf-off conditions. Afterward, three programs were applied to estimate PAI from leaf-on DHP and the woody area index (WAI) from leaf-off DHP. Finally, by subtracting WAI from PAI, true LAI estimates from 37 different algorithms were achieved for evaluation. The performance of these algorithms was compared with pre-defined LAI and PAI values in the virtual forests. The results demonstrated that without correcting for the vegetation clumping effect, Can_Eye, CIMES, and Hemisfer could estimate effective PAI and effective LAI consistent with each other (R2 > 0.8, RMSD < 0.2). After correcting for the vegetation clumping effect, there was a large inconsistency. In general, Can_Eye more accurately estimated true LAI than CIMES and Hemisfer (with R2 = 0.88 > 0.72, 0.49; RMSE = 0.45 < 0.7, 0.94; nRMSE = 15.7% < 24.21%, 32.81%). There was a systematic underestimation of PAI and LAI using Hemisfer. The most accurate algorithm for estimating LAI was identified as the P57 algorithm in Can_Eye which used the 57.5◦ gap fraction inversion combined with the finite-length averaging clumping correction. These results demonstrated the inconsistency of LAI estimates from DHP using different algorithms. It highlights the importance and provides a reference for standardizing the algorithm protocol for in situ forest LAI measurement using DHP.


Introduction
The leaf area index (LAI), defined as one-half of the total green leaf area per unit horizontal ground surface area [1], is a key vegetation structural parameter influencing the process of photosynthesis, transpiration, and rain interception. Because of its importance, LAI has been identified as both an essential climate variable and an essential biodiversity variable [2,3]. In situ LAI measurement plays a vital role in monitoring vegetation dynamics from the ground, as well as for calibrating and validating satellite LAI products. Digital hemispherical photography (DHP) is widely used for in situ LAI measurement. It obtains photographs of the forest vegetation from the ground looking upward through a fisheye lens. By analyzing these photos, the gap fraction can be determined after separating the foliage from the sky, and LAI can be estimated using the gap fraction model. Compared to other techniques such as LAI-2200 and TRAC (tracing radiation and architecture of canopies), DHP has the advantages of lower costs, enhanced visual inspection of canopies, and a permanent archive that can be reprocessed when refined models become available [4]. DHP has been used for vegetation phenology studies [5,6]. DHP collected in the VALERI (Validation of Land European Remote sensing Instruments), and the NEON (National Ecological Observatory Network) was used to evaluate satellite LAI products [7][8][9][10]. LAI estimated from terrestrial [11], airborne [12,13], and spaceborne LiDAR [14] also relied on DHP for validation. Recent studies have pointed out the lack of long-term ground observation and suggested the expansion of existing in situ LAI observation networks as a "high priority" to enhance the quality of satellite-based LAI products [15,16].
There are several steps involved in the processing of DHP to estimate LAI, including the differentiation between sky pixels and canopy pixels, the calculation of gap fraction and gap size, the estimation of clumping index, and the estimation of LAI from gap fraction or gap size distribution inversion [17]. In forests, due to the existence of woody components such as tree trunks and branches, the plant area index (PAI) rather than LAI is estimated from DHP. Effective PAI (PAI eff ) is derived from DHP, assuming that canopy elements are randomly distributed, while true PAI (PAI true ) is derived if the non-randomness of canopy elements is corrected through the estimation of the clumping index.
There are various LAI algorithms implemented in freely available programs to process DHP for PAI/LAI estimation. These algorithms mainly differ in how the PAI eff is estimated, how the orientation and the clumping of leaves are estimated, and how pure canopy segments with zero gap fraction are handled. Some widely used programs include Gap Light Analyzer (GLA) [18], Can_Eye [19], CIMES [4], and Hemisfer [20]. GLA has been continuously used for LAI estimation [5,21,22]. Nevertheless, GLA only provides estimates of PAI eff, as it does not correct for the clumping effect. Can_Eye, developed by the French National Institute of Agricultural Research, has been used extensively in previous studies [23][24][25][26][27]. Hemisfer, developed by the Swiss Federal Institute for the Forest, Snow, and Landscape Research, has been widely used as well [28][29][30]. CIMES is a package of programs encompassing various LAI retrieval methods and is particularly flexible for batch processing multiple DHP images [31][32][33]. Faced with these options, the question that often arises for users is which algorithm of which program produces more accurate LAI estimates. Addressing this problem can provide guidance for standardizing LAI estimation protocols and reducing in situ LAI measurement uncertainty.
Few studies have carried out the accuracy evaluation of different algorithms in commonly used DHP programs when estimating forest LAI. Glatthorn and Beckschäfer compared seven binarization algorithms to classify foliage and sky pixels and found that three algorithms (including Minimum, Edge Detection, and Minimum Histogram) achieved the highest accuracies. This analysis focused only on the gap fraction determination step [34]. Jarčuška et al. (2010) compared the consistency of GLA and WinSCANOPY and found that they produced similar PAI eff estimates [35]. The study by Promis et al. (2011) found similar PAI eff estimates from GLA and Winphot [36]. Similarly, Hall et al. (2017) found that Can_Eye and CIMES produced comparable PAI eff estimates. However, the two programs produced statistically significant different clumping index estimates and, thus, different PAI true estimates [37]. It is worth noting that these studies only evaluated the consistency of results from different algorithms, instead of the accuracy of each algorithm, due to a lack of true reference values. A few studies have used destructive sampling or litterfall collection to acquire LAI reference values for validating the accuracy of DHP in LAI estimation [27,38,39]. However, there is no conclusive evidence concerning the accuracy comparison of various algorithms implemented in these different DHP programs. This calls for an accuracy assessment so that they can be used in the community with confidence.
Since the lack of LAI reference values is the main obstacle for DHP accuracy evaluation, virtual forests offer an alternative platform other than a destructive sampling of all leaves. Virtual forests are a relatively new area still under research. Some researchers define a virtual forest as a computer-based replica of the real forest which is assumed to be of interest for professional and non-professional forest users [40]. Virtual forests can be used for modeling forest growth, predicting forest fire spread, and enabling virtual tourism, as well as calibrating and validating remote sensing data in forest areas [41,42]. There are some previous studies utilizing virtual forests and synthetic DHP to validate the accuracy of the clumping index and PAI true estimates [43,44], leaf angle distribution [45], and slope correction on estimating PAI eff [46]. Nevertheless, these studies used simple geometric primitives to model trees, which differ from real trees, especially in terms of the woody component structure. Recently, highly realistic tree models have been reconstructed from terrestrial LiDAR point clouds and further used to construct virtual forests [47]. Combined with a ray-tracing technique, synthetic DHP is generated for evaluating the retrieval of the clumping index [48]. More recently, Zou et al. (2018) used virtual forests to assess the performance of seven inversion models in estimating the PAI and LAI values from combined leaf-on and leaf-off DHP [49]. Some simulation studies use virtual isolated trees with realistic tree architecture to evaluate the accuracy of leaf area density and LAI estimation for individual trees [50,51]. To the best of our knowledge, there have been no conclusive studies with accuracy evaluations of algorithms implemented in commonly used DHP software programs. Virtual forests provide the potential to solve this problem.
In this study, our research objective is to use virtual broadleaf forests to assess and compare the accuracy of various algorithms implemented in three commonly used DHP software programs in estimating LAI. Both leaf-on and leaf-off virtual forests were created to assess the retrieval of the plant area index (PAI) and leaf area index (LAI). A total of 37 algorithms in three DHP programs, including Can_Eye, CIMES, and Hemisfer, were evaluated. Algorithms and software that do not correct for the clumping effect were not incorporated in the comparison. We aim to provide guidance for users and to identify future directions for the algorithm development of in situ LAI estimation using DHP.

Materials and Methods
The overall experimental design is displayed in Figure 1. Explicit individual tree models (quantitative structure models, QSM) with and without leaves were first constructed by tree reconstruction from point clouds and leaf insertion. Then, the trees were placed randomly to create a series of virtual forest stands with different stem densities and LAI. Synthetic DHP was simulated using a ray-tracing technique, for both leaf-on and leaf-off virtual forests. Afterward, 37 different algorithms utilized in three software packages, including Can_Eye, CIMES, and Hemisfer, were used to process the leaf-on DHP for plant area index (PAI) and leaf-off DHP for woody area index (WAI) estimation. The derived LAI estimates (via subtraction of WAI from PAI) were compared with pre-defined LAI reference values in the virtual forest stands for accuracy evaluation.

Virtual Forests Generation
Realistic tree models were used in this study to construct virtual broadleaf forests. Compared to the simple geometric primitives used in previous studies [43], realistic tree models are used to better simulate the complex structure of forest woody components. For each individual tree, the models for woody components and leaves were generated separately. The process began by reconstructing the 3D models of tree stems and branches with the open-source TreeQSM method as quantitative structure models (QSM) [52]. The TreeQSM method required point cloud data of single leaf-off trees as model inputs. In this study, we used the terrestrial laser scanning point cloud data of European beech (Fagus sylvatica) trees and English oak (Quercus robur L.) trees in leaf-off conditions, which were collected in previous studies [53,54]. A variety of individual trees was constructed, with the heights of 5, 10, 15, 20, 25, and 30 m. A diamond shape consisting of two triangles was used as the base leaf model for all trees. Various leaf sizes were utilized, ranging from 25 to 60 cm 2 . Leaves were inserted to the woody QSM model using a revised non-intersecting leaf insertion algorithm (QSM-FaNNI) [54], so that leaves intersected neither other leaves nor other woody components. Trees of different leaf densities were created. In addition, we revised the original QSM-FaNNI leaf insertion method so that the orientation of all leaves followed pre-defined leaf inclination angle distribution types. As a result, we received a collection of 180 highly realistic tree models, examples of which are shown in Figure 2.

Virtual Forests Generation
Realistic tree models were used in this study to construct virtual broadleaf forests. Compared to the simple geometric primitives used in previous studies [43], realistic tree models are used to better simulate the complex structure of forest woody components. For each individual tree, the models for woody components and leaves were generated separately. The process began by reconstructing the 3D models of tree stems and branches with the open-source TreeQSM method as quantitative structure models (QSM) [52]. The TreeQSM method required point cloud data of single leaf-off trees as model inputs. In this study, we used the terrestrial laser scanning point cloud data of European beech (Fagus sylvatica) trees and English oak (Quercus robur L.) trees in leaf-off conditions, which were collected in previous studies [53,54]. A variety of individual trees was constructed, with the heights of 5, 10, 15, 20, 25, and 30 m. A diamond shape consisting of two triangles was used as the base leaf model for all trees. Various leaf sizes were utilized, ranging from 25 to 60 cm 2 . Leaves were inserted to the woody QSM model using a revised non-intersecting leaf insertion algorithm (QSM-FaNNI) [54], so that leaves intersected neither other leaves nor other woody components. Trees of different leaf densities were created. In addition, we revised the original QSM-FaNNI leaf insertion method so that the orientation of all leaves followed pre-defined leaf inclination angle distribution types. As a result, we received a collection of 180 highly realistic tree models, examples of which are shown in Figure 2.
Afterward, individual tree models were randomly distributed spatially to construct virtual forest stands. A total of 30 scenes were simulated comprising different LAI values and stem distributions as presented in Figure 3. Each forest stand had a size of 120 × 120 m. Each tree was placed with a random rotation around the vertical axis to increase randomness. Rules were defined so that the 3D convex hull of neighboring trees did not intersect with each other, with the implication that a small understory tree could stand beneath a tall tree. A flat topography was assumed for all stands in this study. The detailed stand information is illustrated in Table 1. Afterward, individual tree models were randomly distributed spatially to construct virtual forest stands. A total of 30 scenes were simulated comprising different LAI values and stem distributions as presented in Figure 3. Each forest stand had a size of 120 × 120 m. Each tree was placed with a random rotation around the vertical axis to increase randomness. Rules were defined so that the 3D convex hull of neighboring trees did not intersect with each other, with the implication that a small understory tree could stand beneath a tall tree. A flat topography was assumed for all stands in this study. The detailed stand information is illustrated in Table 1.
The ground-truthed LAI and PAI are termed LAItrue-ref and PAItrue-ref hereafter. They were derived directly from the virtual forest stands, by taking the ratio between "half of the total surface areas of all leaves and all woody components in the forest stand" and "the ground area of the forest stand". To be consistent with the usual height of 1.5 m above ground for most DHP collections, only the leaves and woody components above this height were incorporated in the computation.   (1) Stand PAI true-ref : the reference value of true PAI in the stand (120 m×120 m) from 1.5 m above ground. (2) Stand LAI true-ref : the reference value of true LAI in the stand (120 m×120 m) from 1.5 m above ground.  The ground-truthed LAI and PAI are termed LAI true-ref and PA Itrue-ref hereafter. They were derived directly from the virtual forest stands, by taking the ratio between "half of the total surface areas of all leaves and all woody components in the forest stand" and "the ground area of the forest stand". To be consistent with the usual height of 1.5 m above ground for most DHP collections, only the leaves and woody components above this height were incorporated in the computation.

Synthetic DHP Generation
A ray-tracing tool, the POV-Ray software, was used to generate the synthetic DHP for the virtual forest stands. POV-Ray has been used in previous studies as well [43,45]. For each forest stand, the plot size had a 25 m radius. A regular grid sampling scheme was adopted as suggested by existing literature [17,55]. The specific DHP acquisition locations are displayed in Figure 4. In total, 16 DHP images were acquired for each plot, with the cameras placed 1.5 m above the ground. The image resolution of the synthetic DHP images was 3648 × 3648 pixels. The DHP images were generated for both leaf-on and leaf-off stands. Examples of the synthetic DHP images are presented in Figure 5.

LAI Estimation from DHP
Estimation of LAI from DHP generally relies on the gap fraction inversion. For canopies with randomly distributed leaves, the Poisson model can be used to relate the gap fraction at multiple directions with the amounts and orientations of leaves. For canopies with clumped leaves, models based on the negative binomial probability function were developed [56]. In forests, the effective PAI (PAI eff ) and true PAI (PAI true ) can be related to gap fraction using the equations: where θ is the zenith angle, P(θ) is the gap fraction in the θ direction, G(θ) is the plant projection function value in the θ direction which is determined by the orientation of leaves and woody components [57]. The λ(θ) is the clumping index in the θ direction, which quantifies the degree to which canopy elements deviate from a random distribution. A λ value lower than 1 denotes a clumped distribution. The smaller λ is, the more clumped the canopy is. The determination of G(θ) and λ(θ) are two challenges for PAI true estimation. Using the hinge angle at 57.5 • , G(θ) can be approximated as a constant value of 0.5 regardless of the orientation of plants [55]. It is also possible to use the Miller's integral formula [58] without estimating G(θ): However, it is biased if it is not possible to analyze the DHP in the 0 •~9 0 • range [55,59]. Another common method is to use a one-parameter or two-parameter function to model the leaf angle distribution g(θ) and the G(θ), so as to inverse the G(θ) and PAI from P(θ) observations at multiple directions [60][61][62].
To correct for canopy non-randomness and estimate PAI true [55], three main methods were proposed to estimate the clumping index (λ(θ)), including the finite-length averaging method (LX) [63], the gap size distribution method (CC) [64,65], and the combination of LX and CC (CLX) method [66]. The finite-length averaging (LX) method was proposed by Lang and Xiang in 1986 using the following equation: where P(θ) is the canopy mean gap fraction of all segments, and ln P(θ) is the logarithmic mean from gap fractions of all segments [63]. However, two assumptions underlie this method, i.e., the foliage within the finite length segment is random, and the segment contains gaps. For segments completely obstructed by leaves, both Can_Eye and CIMES adopted a saturated LAI (L sat ) value to address this problem. In this experiment, L sat was set as the default value of 10 following the manuals of Can_Eye and CIMES. The gap size distribution method was proposed by Chen and Cihlar in 1995 and corrected by Leblanc in 2002, using the following equation: where F m (0, θ) is the accumulated canopy gap fraction, F mr (0, θ) is the reduced gap size accumulated fraction after removal of large non-random gaps from the measured gap size accumulation curve F m (λ) until the pattern of F mr (0, θ) resembles that of an equivalent canopy with a random spatial distribution of foliage [64,65]. In 2005, Leblanc proposed to combine the LX and the CC methods to address the potentially violated assumption of a random distribution of canopy elements at the segment scale associated with the LX method by using: where P k (θ) is the gap fraction of the segment k, and λ CCk (θ) is the clumping index of the segment k using the CC method. Another method to estimate PAI true consists of averaging values of local PAI eff over azimuthal angular intervals (WT method) [67]. Since forests contain many woody components other than leaves, a woody effect correction is necessary to convert the PAI true to LAI true . In this study, we used leaf-on and leaf-off DHP to estimate PAI and woody area index (WAI), respectively [68]. Then, the LAI was estimated using: where PAI true is the estimate of PAI after considering canopy non-randomness in leaf-on conditions, WAI true is the estimate of WAI (PAI in leaf-off conditions) after considering canopy non-randomness, and LAI true is the final estimate of LAI. Specifically, all DHP images for both leaf-on and leaf-off virtual forests in this study were firstly processed with the automatic two-corner method to separate the sky from canopy pixels [69]. The two-corner method proved to be stable and more accurate than other methods based on previous studies [70]. After the image classification, the binary DHP images were imported into all DHP software, including Can_Eye, CIMES, and Hemisfer, for LAI estimation. In all software programs, the DHP was broken into multiple sub-sectors, with a 2.5 • zenith angular resolution and a 10 • azimuth angular resolution. The sub-sectors with zenith angles above 60 • were removed from subsequent analysis due to a high portion of mixed pixels. All software are regularly improved and updated. The Can_Eye software was used is the version 6.495. The Hemisfer software was used is the version 3. For CIMES, we used the version (1982-2020).
In total, there were 37 algorithms from Can_Eye, CIMES, and Hemisfer producing PAI true and LAI true estimates (4 from Can_Eye, 9 from CIMES, and 24 from Hemisfer). The differences among these algorithms were mainly in how the PAI eff was estimated, how the orientation of leaves and the clumping index were estimated, and how pure canopy segments with zero gap fraction were handled. Summarized descriptions of the 37 algorithms are presented in Tables 2-4 for more details.
In general, Can_Eye offered estimates of PAIeff and PAItrue using either a single direction (57.5 • ) gap fraction or multidirectional (0 • −60 • ) gap fraction with different inversion methods. The clumping index λ in Can_Eye was estimated based on the LX method [63]. Dissimilar to Can_Eye, CIMES was able to estimate λ not only using the LX but also using the CC and CLX methods [4]. The Hemisfer program implemented the LX and CC methods for clumping correction, as well as the non-linearity gap fraction correction. It is worth noting that even when using the same basic algorithm, the detailed procedures and inversion schemes may be inconsistent between different software. It is suggested to refer to the user manual of each software program for a more detailed description of each algorithm.

Algorithm Abbreviation
Basic Principle References P57 1. Use of the gap fraction at 57.5 • (55 •~6 0 • ) 2. The G(θ) was approximated as 0.5 regardless of g(θ) types 3. Clumping correction was based on the LX method 4. L sat was set as 10 for pure segments with no gaps [19,63] v5.1 1. Use of the gap fraction at 0 •~6 0 • 2. The g(θ) which determined G(θ) was modeled by the ALA using the ellipsoidal distribution 3. PAI and ALA were inversed using a look-up table scheme, with the cost function constrained by a term of ALA (the retrieved ALA value must be close to 60 • ± 30 • ) 4. Clumping correction was based on the LX method 5. L sat was set as 10 for pure segments with no gaps [19,60,63] v6.1 1. Use of the gap fraction at 0 •~6 0 • 2. The g(θ) was modeled by the ALA using the ellipsoidal distribution 3. PAI and ALA were inversed using a lookup table scheme, with the cost function constrained by a term of PAI 57 (the retrieved PAI value that must be close to the one retrieved from the annulus at 57.5 • ) 4. Clumping correction was based on the LX method 5. L sat was set as 10 for pure segments with no gaps [19,60,63] Miller 1. Use of the gap fraction at 0 •~6 0 • 2. Use of Miller's formula to estimate PAI eff 3. Clumping correction based on the LX method 4. L sat was set as 10 for pure segments with no gaps [19,58]

Statistical Analysis
In terms of the PAI eff and LAI eff , we calculated the consistency among the three programs (Can_Eye, CIMES, and Hemisfer), using the coefficient of determination (R 2 ), and the root mean square difference (RMSD). Higher values of R 2 and lower values of RMSD indicated greater consistency and robustness.
As for PAI true and LAI true , the values calculated from the virtual forests as described in Section 2.1 were used as the true reference values. Because the three programs offered 37 estimates of PAI and LAI from different algorithms, we first identified the most accurate results within each software program. In addition, the most accurate results were subsequently used for inter-software comparison, in terms of R 2 , RMSE, and normalized RMSE (nRMSE). Higher values of R 2 , lower values of RMSE, and lower values of nRMSE indicate higher accuracy.

PAI eff and LAI eff Estimation Results
Without correcting for the clumping effect caused by vegetation non-randomness, the estimates of PAI (PAI eff-est ) were on average 55.8% of the PAI true-ref values, while the estimates of LAI (LAI eff-est ) were on average 51.22% of LAI true-ref values. The comparison of the PAI eff-est and LAI eff-est among the three programs is shown in Figure 6. The PAI eff and LAI eff estimates from CIMES were slightly higher than those from Can_Eye were (as shown in Figures 6a,d. Compared to CIMES, the PAI eff and LAI eff estimates from Hemisfer were closer to Can_Eye (RMSD = 0.11 < 0.19 for PAI eff-est , and RMSD = 0.09 < 0.14 for LAI eff-est , as shown in Figures 6b,e. In general, the results of effective PAI and effective LAI from the three programs were consistent (R 2 ≥ 0.8, RMSD ≤ 0.19).

Comparison of PAItrue Estimation Accuracy
Using the PAIeff-est from Can_Eye divided by the plot PAItrue-ref, we obtained the clumping index values (λ) of each forest plot, as shown in Table 5. This quantified the level of the clumping effect in each stand and assisted in the evaluation of clumping correction methods using different algorithms.

Comparison of PAI true Estimation Accuracy
Using the PAI eff-est from Can_Eye divided by the plot PAI true-ref , we obtained the clumping index values (λ) of each forest plot, as shown in Table 5. This quantified the level of the clumping effect in each stand and assisted in the evaluation of clumping correction methods using different algorithms.
All results of PAI true-est using different algorithms in Can_Eye, CIMES, and Hemisfer are presented in Figures 7-9, respectively. The symbol of each plot in each figure was colored based on the clumping index (λ) value according to the results in Table 5, and the size was adjusted according to the average leaf inclination angle of each plot in Table 1.  Figure 7). The most accurate algorithm was the P57 algorithm, which used the gap fraction (at 57.5°) inversion combined with the LX clumping correction method (R 2 = 0.86, RMSE = 0.49, nRMSE = 13.64%). The least accurate algorithm in Can_Eye was Miller's formula using the gap fraction at 0°~60° combined with the LX clumping correction method (RMSE = 1.68, nRMSE = 46.24%, see Figure 7). In terms of CIMES, the nine algorithms produced quite different PAItrue-est values, with nRMSE ranging in (19.3%, 54.37%) (see Figure 8). The most accurate algorithm of CIMES was from the multiple direction gap fraction inversion at 0°~60° using the Campbell approach combined with the LX clumping correction (CAM_LX algorithm, R 2 = 0.73, RMSE = 0.7, nRMSE = 19.3%), while the least accurate algorithm was the Mil-ler_CC57 method, with the nRMSE reaching 54.37% (see Figure 8). When comparing PAI true-est with PAI true-ref , the four algorithms in Can_Eye produced different PAI true-est values, with nRMSE ranging in (13.64%, 46.24%) (see Figure 7). The most accurate algorithm was the P57 algorithm, which used the gap fraction (at 57.5 • ) inversion combined with the LX clumping correction method (R 2 = 0.86, RMSE = 0.49, nRMSE = 13.64%). The least accurate algorithm in Can_Eye was Miller's formula using the gap fraction at 0 •~6 0 • combined with the LX clumping correction method (RMSE = 1.68, nRMSE = 46.24%, see Figure 7).
In terms of CIMES, the nine algorithms produced quite different PAI true-est values, with nRMSE ranging in (19.3%, 54.37%) (see Figure 8). The most accurate algorithm of CIMES was from the multiple direction gap fraction inversion at 0 •~6 0 • using the Campbell approach combined with the LX clumping correction (CAM_LX algorithm, R 2 = 0.73, RMSE = 0.7, nRMSE = 19.3%), while the least accurate algorithm was the Miller_CC57 method, with the nRMSE reaching 54.37% (see Figure 8).
An inter-comparison among all 37 algorithms revealed that the most accurate algorithm to estimate PAI true was the P57 method in Can_Eye, which used the gap fraction (at 57.5 • ) inversion combined with the LX clumping correction. Of note, there was a strong systematic underestimation of PAI true by Hemisfer, either for canopies with high or low clumping. The PAI true estimates from Hemisfer reached saturation at PAI values around three (Figure 9j). In comparison, the P57 method in Can_Eye could accurately correct the clumping effect until the vegetation reached high clumping (when λ < 0.45); then, it began to underestimate PAI true in forests with PAI above 4.5 (Figure 7d). There was neither a systematic underestimation nor overestimation of PAI true using the CAM_LX algorithm from CIMES (Figure 8a). In Figures 7d, 8a and 9j, there was no clear spatial distribution pattern of symbol sizes, implying that the average leaf inclination angle of each stand had little effect on the estimation of PAI true . Regarding Hemisfer, the 24 algorithms produced different PAItrue-est values, with nRMSE ranging in (30.46%, 43.34%) (see Figure 9). The most accurate PAItrue-est result from Hemisfer was obtained with the LX_Miller method (R 2 = 0.32, RMSE = 1.11, nRMSE = 30.46%), while the least accurate algorithm was the CC_2000 method, with the nRMSE reaching 43.34% (see Figure 9).
An inter-comparison among all 37 algorithms revealed that the most accurate algorithm to estimate PAItrue was the P57 method in Can_Eye, which used the gap fraction (at 57.5°) inversion combined with the LX clumping correction. Of note, there was a strong systematic underestimation of PAItrue by Hemisfer, either for canopies with high or low clumping. The PAItrue estimates from Hemisfer reached saturation at PAI values around three (Figure 9j). In comparison, the P57 method in Can_Eye could accurately correct the clumping effect until the vegetation reached high clumping (when λ < 0.45); then, it began to underestimate PAItrue in forests with PAI above 4.5 (Figure 7d). There was neither a systematic underestimation nor overestimation of PAItrue using the

Comparison of LAI true Estimation Accuracy
In terms of LAI true-est , similar to the case of PAI true-est , different algorithms produced quite different LAI true-est results (nRMSE ranged in (15.7%, 53.24%) for Can_Eye, (24.21%, 70.64%) for CIMES, and (32.81%, 49.49%) for Hemisfer). Within each software program, the most accurate algorithm to estimate LAI true was the same as PAI true . For more details, the reader can refer to Figures S1-S3 in the Supplementary. In the following, only the most accurate algorithm in each software program was listed for an intercomparison.
The P57 method in Can_Eye which used the gap fraction (at 57.5 • ) inversion combined with the LX clumping correction was revealed as the most accurate algorithm to estimate LAI true compared to CIMES and Hemisfer (R 2 = 0.88 > 0.72, 0.49; RMSE = 0.45 < 0.7, 0.94; nRMSE = 15.7% < 24.21%, 32.81%, see Figure 10). There was a more severe underestimation of LAI true from the LX_Miller algorithm in Hemisfer than from CIMES and Can_Eye, even in forests with a moderate amount of leaves at the LAI value around 2.5 ( Figure 10c). In comparison, the P57 method in Can_Eye started to underestimate LAI in forests with dense leaves, with an LAI value around four (Figure 10a). There was neither a systematic underestimation nor overestimation of LAI true from the CAM_LX algorithm in CIMES (Figure 10b).

Comparison of LAItrue Estimation Accuracy
In terms of LAItrue-est, similar to the case of PAItrue-est, different algorithms produced quite different LAItrue-est results (nRMSE ranged in (15.7%, 53.24%) for Can_Eye, (24.21%, 70.64%) for CIMES, and (32.81%, 49.49%) for Hemisfer). Within each software program, the most accurate algorithm to estimate LAItrue was the same as PAItrue. For more details, the reader can refer to Figures S1, S2, and S3 in the Supplementary. In the following, only the most accurate algorithm in each software program was listed for an intercomparison.
The P57 method in Can_Eye which used the gap fraction (at 57.5°) inversion combined with the LX clumping correction was revealed as the most accurate algorithm to estimate LAItrue compared to CIMES and Hemisfer (R 2 = 0.88 > 0.72, 0.49; RMSE = 0.45 < 0.7, 0.94; nRMSE = 15.7% < 24.21%, 32.81%, see Figure 10). There was a more severe underestimation of LAItrue from the LX_Miller algorithm in Hemisfer than from CIMES and Can_Eye, even in forests with a moderate amount of leaves at the LAI value around 2.5 ( Figure 10c). In comparison, the P57 method in Can_Eye started to underestimate LAI in forests with dense leaves, with an LAI value around four (Figure 10a). There was neither a systematic underestimation nor overestimation of LAItrue from the CAM_LX algorithm in CIMES (Figure 10b).

Discussion
The results of this study demonstrated that three commonly used programs for DHP (including Can_Eye, CIMES, and Hemisfer) could estimate consistent effective PAI and effective LAI, though these programs were inconsistent when estimating true PAI and true LAI using the same DHP images (Figures 6-10). Our results are in agreement with a previous study by Hall et al (2017), who also found similar PAIeff estimates, but significantly different estimates of PAItrue between Can_Eye and CIMES [37]. However, this previous study by Hall et al. only reported the inconsistency between Can_Eye and CI-MES. The contribution of our research is that we proved that Can_Eye estimates more accurate PAItrue and LAItrue values than CIMES and Hemisfer. In addition, we identified the most accurate algorithm among all 37 algorithms as the P57 algorithm in Can_Eye.

Discussion
The results of this study demonstrated that three commonly used programs for DHP (including Can_Eye, CIMES, and Hemisfer) could estimate consistent effective PAI and effective LAI, though these programs were inconsistent when estimating true PAI and true LAI using the same DHP images (Figures 6-10). Our results are in agreement with a previous study by Hall et al (2017), who also found similar PAI eff estimates, but significantly different estimates of PAI true between Can_Eye and CIMES [37]. However, this previous study by Hall et al. only reported the inconsistency between Can_Eye and CIMES. The contribution of our research is that we proved that Can_Eye estimates more accurate PAI true and LAI true values than CIMES and Hemisfer. In addition, we identified the most accurate algorithm among all 37 algorithms as the P57 algorithm in Can_Eye. This can guide future users during the algorithm and software selection step when measuring in situ LAI from DHP. In addition, the use of virtual forests was not subject to inherent field measurement errors compared to litterfall collection or destructive sampling. It, thus, provides strong confidence in the results of the algorithm assessment.
The results of this study indicate that the major difference in estimating PAI true and LAI true among different software programs (i.e., Can_Eye, CIMES, and Hemisfer) was due to different estimates of the clumping index. Many researchers have pointed out that the estimates of the clumping index remained a large source of uncertainty for LAI estimation [44,75]. There was a systematic underestimation of PAI true and LAI true using Hemisfer, even in forests with a moderate amount of leaves and low vegetation clumping (as shown in Figure 9g-l and Figure S3g-l). This may have been caused by the pure canopy segments with no gaps. In such cases, the underlying assumption "the segment contains gaps" in the LX method did not hold. Both the Can_Eye and CIMES software utilized a saturated LAI (L sat ) value to address this problem. In this experiment, L sat was set as 10 for both Can_Eye and CIMES. Dissimilar to Can_Eye and CIMES, Hemisfer did not provide this L sat solution and offered no details about how pure canopy segments were handled. This may be the reason for the systematic underestimation of PAI true and LAI true in Hemisfer.
The better performance of the P57 algorithm in Can_Eye was perhaps because it used the 57.5 • angle when G(θ) was approximated as 0.5, while other methods such as the v5.1 algorithm and the v6.1 algorithm in Can_Eye and the CAM_LX algorithm in CIMES used multiple directions (0 •~6 0 • ) and had to inverse both the average leaf inclination angle (ALA) as well as PAI or LAI, which could add more uncertainty and errors (as shown in Figures 7b-c, 8a, S1b-c and S2a). Even using the same 0 •~6 0 • direction and the same LX clumping correction method, CIMES generated a less accurate PAI true and less accurate LAI true than Can_Eye (i.e., comparing Figure 8a with Figure 7b,c and comparing Figure  S1a with Figure S2b,c in the Supplementary). This may be caused by the different inversion schemes implemented in Can_Eye and CIMES. Can_Eye introduces a regularization term in the cost function, which imposes constraints to improving the PAI estimates (see Equations (8) and (9)), while CIMES is not similarly constrained [4,19,59]. The higher accuracy of the v6.1 method than the v5.1 method in Can_Eye found in our study (R 2 = 0.86 > 0.81, RMSE = 0.5 < 0.58, nRMSE = 13.66% < 15.93%, as shown in Figure 7) corroborates the prediction of higher efficiency of the PAI regularization term when compared with the ALA regularization term [19].
v5.1 method : In terms of the CC clumping correction method [64,65], there was a clear underestimation of PAI true and LAI true using CIMES and Hemisfer (as shown in Figure 8f,g, Figure 9a-f,m-r, Figures S2f,g and S3a-f,m-r in the Supplementary). This is in agreement with previous studies, which also found that the CC method underestimated the clumping effect (i.e., overestimation of λ [32,76]. Previous studies [32,43,66] suggested the CLX method worked better than the LX and CC method because it addressed the issues of "clumping within the segment". Currently, only CIMES implemented the CLX clumping correction method. Therefore, this study cannot present a comparison among different software implementations of the CLX methods. However, in our study, we found that the CLX method using multiple directions (0 •~6 0 • ) in CIMES generally failed to estimate the WAI using leaf-off DHP (see Figure S2i in the Supplementary). Further inspection revealed that the CLX estimate using the LAIMIL method in CIMES equaled 0. Currently, we cannot offer a clear explanation about this issue. It is suspected that there may have been some errors in the implementation of the CLX correction in the LAIMIL method. Nevertheless, it is suggested to be cautious using the CLX method in CIMES before its refinement and perhaps to incorporate the CLX method in Can_Eye, which may result in more accurate estimates of PAI true and LAI true .
The woody component was another obstacle for calculating the true LAI. One method to correct the woody effect was to multiply the estimate of PAI by a woody-to-total area ratio (α) to achieve WAI. An empirical α value could be estimated based on destructive sampling for a certain type of tree. Such manually intensive methods are dependent on species composition and growth form, and may not capture the spatial variation of α across stands [77]. The use of a near-infrared (NIR) channel camera to differentiate between woody and foliar components has been suggested but has not been applied widely [78]. Some studies combined the use of DHP acquired in leaf-on and leaf-off conditions for broadleaf forests, to estimate LAI by subtracting leaf-off PAI from leaf-on PAI [68], though the results were not validated. The study by  proved that subtraction of WAI eff from PAI eff correlated well with LAI eff when using a virtual forest [79]. The results of our study further illustrated that the subtraction of WAI true from PAI true correlated well with pre-defined LAI true values ( Figure 10). This demonstrated the effectiveness of using leaf-on and leaf-off DHP to correct the woody effect and estimate LAI for broadleaf forests. This can also provide a reference for dynamic LAI true monitoring studies using repeated DHP for broader environmental applications [80].
Continuous and freely available satellite LAI products have boosted global Earth system modeling studies [15]. However, an inconsistency was found in different satellite LAI products such as GLASS, LAI3g, TCDR, and GLOBMAP [81]. Many more in situ LAI measurements are in great need to enhance the accuracy of satellite LAI products [1]. Previous studies have highlighted the uncertainty in in situ LAI measurements using different techniques, including DHP, terrestrial LiDAR, and LAI-2200 [70]. Our study highlighted the inconsistency of LAI estimates from DHP using different algorithms from various software programs. It proves that care must be taken, and analysts should carefully follow a methodological protocol while processing the DHP for LAI estimates. The finding of this study can guide the expansion of ground LAI observation networks, which is vital to enhance the quality of satellite products [16]. Further studies must be implemented to consolidate the accuracy of in situ LAI measurements, which is the basis for calibrating and validating satellite LAI products for downstream applications.
Virtual forests can play an important role in refining in situ LAI estimation methods, since high-fidelity, ground-reference LAI, such as that from destructive measurements, is rarely available at stand level. To acquire true values of LAI and PAI in real forests, we would need to destructively sample all leaves and woody components and measure their surface area in forests. This is too expensive and impractical. Although the litterfall collection method can be used to estimate true LAI, it is prone to uncertainty caused by litter trap frequency and collection time intervals [82]. However, in the set-up of virtual forests, we can control and calculate all aspects ranging from each leaf to each branch of the forests. The virtual forest created in this research had PAI values ranging from 1.43 to 6.38, which is consistent with a previous study that shows a range of PAI values from 1.1 to 8.0 in temperate deciduous broadleaf forests [83]. In addition, the clumping index ranging from 0.38 to 0.78 in this study is also similar to the range of 0.41 to 0.92 for clumping indices of forests derived from MODIS data in previous studies [84,85]. This strengthens the reliability of the virtual forest generation framework employed in this research. Further studies may adopt the actual position of individual trees from the forest inventory instead of random displacements of trees. A recent approach named "Canopy Constructor", which combines field inventories and airborne LiDAR scans to create virtual forests that best fit the data, is a promising method [86]. Virtual forests can be used in the algorithm development stage as well as indoor evaluation at low cost. It offers the possibility for DHP re-collection, which is useful when a new LAI collection protocol is designed. Although field measurement is still indispensable, virtual forests can serve as a supplement for validating the theoretical error of each factor inherent in LAI measurement and retrieval [87].
Several limitations can be addressed in future relevant research. Firstly, the current study assumed a relatively small forest size due to the limitation of computing resources. Secondly, for simplicity, our study assumed a flat terrain and pure broadleaf stands in all virtual forests. Given that much progress has been achieved in reconstructing highly realistic coniferous trees from terrestrial LiDAR point clouds [88,89], future studies can investigate the accuracy of in situ LAI measurement in virtual coniferous or mixed forests and mountainous forests.

Conclusions
DHP is a widely used in situ LAI technique, which was used to collect LAI ground references for improving satellite LAI products. Given the array of various algorithms and software programs that are available for processing DHP for LAI estimation, an accuracy evaluation among these methods has become urgent. There have been few studies comparing the accuracy of these algorithms in forests due to the difficulty in achieving a high-fidelity reference LAI value. This study used virtual forests to address this issue. Various broadleaf forests covering a range of tree heights, stem densities, and LAI values were constructed. Combined with ray-tracing, a synthetic DHP was simulated, both for leaf-on and leaf-off conditions. A total of 37 algorithms from three programs, including Can_Eye, CIMES, and Hemisfer, were used to estimate the plant area index (PAI) from leaf-on DHP, woody area index (WAI) from leaf-off DHP, and LAI by subtracting WAI from PAI. Other software such as the Gap Light Analyzer and HemiView were not investigated in this study since they do not correct for the clumping effect. Overall, the results demonstrated that effective PAI and LAI results could be consistently retrieved from the three programs (R 2 > 0.8, RMSD < 0.2). However, after correcting for the clumping effect caused by a non-random vegetation distribution, a large inconsistency occurred between the estimates of true PAI and true LAI from different algorithms. The results demonstrated that Can_Eye more accurately estimated true LAI than CIMES and Hemisfer (with R 2 = 0.88 > 0.72, 0.49; RMSE = 0.45 < 0.7, 0.94; nRMSE = 15.7% < 24.21%, 32.81%). It was also found that the P57 algorithm, which used the 57.5 • gap fraction inversion combined with the finite-length averaging (LX) clumping correction method implemented in Can_Eye, was the most accurate. This study highlighted the in situ LAI measurement uncertainty due to the choice of the LAI algorithm and demonstrated the effectiveness of using leaf-on and leaf-off to correct for the woody component's effect on LAI estimation. From the results of this study, we recommend using the 57.5 • gap fraction inversion combined with the finite-length averaging (LX) clumping correction method implemented in Can_Eye for LAI estimation in broadleaf forests. Further studies exploring coniferous and mixed forests are suggested to consolidate the protocol of in situ LAI measurements to reduce in situ LAI uncertainty.
Supplementary Materials: The following supplementary materials are in the Supplementary available online at https://www.mdpi.com/article/10.3390/rs13163325/s1, Figure S1: LAI true results using different algorithms from Can_Eye, Figure S2: LAI true results using different algorithms from CIMES, Figure S3: LAI true results using different algorithms from Hemisfer.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.