Extracting Leaf Area Index by Sunlit Foliage Component from Downward-Looking Digital Photography under Clear-Sky Conditions

The development of near-surface remote sensing requires the accurate extraction of leaf area index (LAI) from networked digital cameras under all illumination conditions. The widely used directional gap fraction model is more suitable for overcast conditions due to the difficulty to discriminate the shaded foliage from the shadowed parts of images acquired on sunny days. In this study, a new LAI extraction method by the sunlit foliage component from downward-looking digital photography under clear-sky conditions is proposed. In this method, the sunlit foliage component was extracted by an automated image classification algorithm named LAB2, the clumping index was estimated by a path length distribution-based method, the LAD and G function were quantified by leveled digital images and, eventually, the LAI was obtained by introducing a geometric-optical (GO) model which can quantify the sunlit foliage proportion. The proposed method was OPEN ACCESS Remote Sens. 2015, 7 13411 evaluated at the YJP site, Canada, by the 3D realistic structural scene constructed based on the field measurements. Results suggest that the LAB2 algorithm makes it possible for the automated image processing and the accurate sunlit foliage extraction with the minimum overall accuracy of 91.4%. The widely-used finite-length method tends to underestimate the clumping index, while the path length distribution-based method can reduce the relative error (RE) from 7.8% to 6.6%. Using the directional gap fraction model under sunny conditions can lead to an underestimation of LAI by (1.61; 55.9%), which was significantly outside the accuracy requirement (0.5; 20%) by the Global Climate Observation System (GCOS). The proposed LAI extraction method has an RMSE of 0.35 and an RE of 11.4% under sunny conditions, which can meet the accuracy requirement of the GCOS. This method relaxes the required diffuse illumination conditions for the digital photography, and can be applied to extract LAI from downward-looking webcam images, which is expected for the regional to continental scale monitoring of vegetation dynamics and validation of satellite remote sensing products.


Introduction
Leaf area index (LAI) is a key biophysical variable to characterize vegetation canopy structure and functioning in most ecosystem productivity and land surface process models [1,2].LAI is defined as half the total foliage area per unit ground surface area [3].Satellite remote sensing provides the unique way to obtain LAI in long-term time series and at the global scale [2,4].However, the accuracy of remotely sensed LAI can be affected by the land surface heterogeneities, the impact of clouds and aerosols in the atmospheric correction, the uncertainties from the forward model used to create the look-up tables, and the saturation of optical signals over dense canopies when the lower layer is obscured by the upper layer [5][6][7][8].Thus, it is necessary to validate the remotely sensed LAI products with ground-based LAI measurements for product utilization and algorithm improvement.In general, LAI can be measured through direct and indirect methods in the field campaign [9].Direct LAI measurements include destructive sampling and non-harvest litter traps for deciduous forests, which are the most accurate, but are extremely labor-intensive and time-consuming [9][10][11].Indirect methods using optical radiometric or imaging sensors, e.g., LAI-2000 Plant Canopy Analyzer (PCA), Digital Hemispherical Photography (DHP), and Tracing Radiation and Architecture of Canopies (TRAC), are based on the gap fraction or gap size distribution analysis [9,[12][13][14].Indirect LAI measurements are widely used in the field campaign due to the high efficiency, but the temporal revisit frequency and the spatial coverage are limited by the manpower [4,15].
Recently, near-surface remote sensing using networked digital cameras or radiometric sensors provides a low-cost way to continuously monitor the vegetation dynamics at high temporal frequency (several measurements per day) over the regional to continental scale [15][16][17][18][19][20].For example, the PhenoCam dataset has collected time series of red-green-blue (RGB) camera images for more than 200 forest sites across North America [16,20].The downward inclined cameras are installed at the top of instrument towers, such as the flux towers for eddy covariance measurements, which are several or tens of meters above the canopy [17].The cameras are set at automated or fixed exposure, and webcam images are widely used to extract color indices in long-term time-series, such as the excess green (ExG) and the green chromatic coordinate (gcc) for vegetation phonological monitoring [16,17,20].Recently much attention has been paid to use LAI time-series for tracking vegetation dynamic changes, because color indices or spectral indices, which are derived from camera RGB channel digital numbers (DNs) or reflectances, can be affected by leaf optical properties (e.g., chlorophyll a+b content and water amount), canopy structure (e.g., LAI, leaf angle distribution, and clumping index), soil background, sun-target-sensor geometry, and diffuse irradiance ratio [21][22][23].Additionally, compared with vegetation indices, the quantitative biophysical meaning of LAI is clear, and LAI will not get saturated on the mature phase at least by definition [24][25][26].However, little research has been conducted to extract LAI from such downward-looking webcam images, which is expected for the regional to continental scale monitoring of vegetation dynamics and validation of satellite remote sensing products.
Current indirect LAI measurements by upward-pointing or downward-looking imaging sensors, such as DHP or non-fisheye digital photography, require stable diffuse illumination conditions, including uniform overcast, just before sunrise or after sunset [15,27].Downward-looking cameras have already been used to extract LAI over short vegetation, such as agricultural crops, but are not often used for tall forest canopies until the tower-based webcam provides the possibility to extract forest LAI by near-surface remote sensing [27][28][29][30].The series of acquired webcam images can be not only under overcast conditions, but also under clear-sky conditions with direct sunlight [16,17].According to the theory of the geometric-optical (GO) model, the difference between the two types of illumination conditions for the images is that on overcast days there are two scene components: foliage and background, while on sunny days there are four scene components: sunlit foliage, sunlit background, shaded foliage, and shaded background [31][32][33].The directional gap fraction model has been widely used on overcast days, while till now little research is dedicated to the extraction of LAI on sunny days [29,30].
The main challenge for the images acquired on sunny days is that it is difficult to discriminate the shaded foliage and shaded background from the shadows, where the signal variations in the three color channels are small [28][29][30].Misclassification of shaded foliage as background on sunny days will lead to an overestimation of the gap fraction extraction [29], which may eventually result in the LAI estimation be outside the relative accuracy (20%) and absolute accuracy (0.5) requirement by the Global Climate Observation System (GCOS, http://www.wmo.int/pages/prog/gcos).To avoid the potential impact of such misclassification on sunny days, the goal of this study is to extract LAI by the area ratio of the sunlit foliage component using the GO model, instead of using the directional gap fraction model for images acquired under overcast conditions.The sunlit foliage component extraction, clumping index estimation and the LAI retrievals by different methods will be presented in detail in this study.The accurate extraction of LAI from downward-looking digital photography under clear-sky conditions will make it possible to monitor the LAI dynamics from low-cost webcams under all illumination conditions.

Theory
According to the GO model, the digital image of the scene within the field of view can be divided into four components: sunlit foliage, sunlit background, shaded foliage, and shaded background [31][32][33].PT, PG, ZT, and ZG are the area ratio of the four components, respectively, and thus PT + PG + ZT + ZG = 1 holds.The area ratio of the sunlit foliage PT for discontinuous canopies can be expressed by [33]: where μ  = θ  , μ  = θ  , and θ  and θ  are the solar/view zenith angle, respectively.Ω  and Ω  are the clumping index in the solar and view directions [12], and   and   are the mean projection of unit foliage area on a plane perpendicular to the solar and view directions [34].The hot spot function w can be derived by: where H is the canopy height, and   is the characteristic linear dimension of the foliage in the Kuusk Hot-spot model [35].For the spherical orientation of leaves, the foliage dimension   =    2 /16, and for horizontal leaves   =   /4, where   is the mean foliage diameter [36].δ = √ From Equations ( 1)-( 3), LAI can be derived by the area ratio of the sunlit foliage component PT as: where the variables have the meaning as in Equation ( 1)-( 3).Thus, the following procedures are proposed to extract LAI from downward-looking digital images acquired on sunny days (Figure 1): 1. Estimate the area ratio of the sunlit foliage component PT from digital images by an automated image classification algorithm.2. Extract the clumping index Ω from digital images by a path length distribution-based method [37].3. Characterize the leaf angle distribution (LAD) and calculate the leaf projection function (G). 4. Acquire the canopy height H, the foliage diameter   and the solar/view geometric information by field measurements.5. Derive LAI by Equation (4) with the above variables estimated.
The procedures (1-2) will be described in detail at the following Sections 2.2-2.3.

Sunlit Foliage Component
An automated image classification algorithm called LAB2 is adopted in this study to extract sunlit foliage from images acquired on sunny days [38].The LAB2 method was originally proposed to estimate the foliage cover from digital nadir-view images; and has shown to be the most accurate method for foliage cover > 10% when compared with other five classification methods [38].
Firstly, the green leaf algorithm (GLA) for each pixel is calculated as [39] GLA = 2 −  −  2 +  + where R, G, and B are the DNs of the RGB channels of the images.The values of GLA range between −1 and +1, and the positive values indicate green vegetation because they have higher levels in the green channel than in the red and blue channels.Secondly, the digital image is converted from the RGB color space to the CIE L*a*b* color space, where L* is the luminance component, a* represents color on the green-magenta axis, and b* represents color on the blue-yellow axis.The CIE L*a*b* color space separates the luminance channel L* from two chromaticity channels a* and b*, which makes the correlations between channels minimal, and reduces the potential impact of illumination changes on image processing, compared with the RGB space.Negative a* indicates green, while positive a* indicates red, which has been implemented in previous studies to quantify the greenness of foliage [38,40].
Thirdly, automatically identify the background and sunlit foliage training sets by the DNs of the RGB images.Pixels with GLA ≤ 0 are classified as definite background training sets, while pixels satisfying the criterion (G > R) & (G > B) & (G > 25) are identified as sunlit foliage training sets, according to the suggestions by [38] to detect the background and foreground training sets.Mean values of GLA, a*, and b* are calculated for the background and sunlit foliage training sets.
Finally, classify each pixel by the Euclidean distance from the sunlit foliage and background means of GLA, a*, and b* using the minimum-distance-to-means classifier.The area ratio of the sunlit foliage PT can be calculated by NT/N, where N is the total number of pixels for the digital image, and NT is the number of pixels classified as the sunlit foliage.

Clumping Index
The clumping index (CI), which quantifies the degree of foliage nonrandom distribution in space, is defined as [12]: where   is the true LAI of the scene, and   (θ  , φ  ) is the effective LAI that can be derived from the directional gap fraction model: where (θ  , φ  ) is the directional gap fraction at the view zenith angle θ  and view azimuth angle φ  , (θ  , φ  ) is the foliage projection function as in Equation (1).(θ  , φ  ) can be extracted from images acquired on adjacent overcast days by the LAB2 image classification algorithm in Section 2.2.However, in practice, the true LAI is usually unknown, which means the clumping index cannot be directly derived by definition as in Equation (6).The total clumping index Ω can be divided into two components where Ω  is the element clumping index quantifying the degree of foliage clumping at scales larger than the shoot, and γ  is the needle-to-shoot area ratio, which is usually assumed to be 1 for broadleaf forests [14].In general, the element clumping index Ω  can be estimated by the finite-length logarithmic gap averaging (LX) method of [41] or gap size distribution-based (CC) method of [42] from sampled transects in the image.In this study, a path length distribution-based method is adopted to estimate the clumping index since it has shown to improve the accuracy on the basis of the widely used LX and CC methods over heterogeneous canopies [37].The estimated true LAI in Equation ( 6) by the path length-based method can be expressed as [37]: where  is the foliage area volume density,   is the maximum path length along the transect, θ  is the view zenith angle,  is the relative path length, and   () is the path length distribution function inversed from measured gaps in the sliding windows.Similar to the segment length in the LX method, the size of the sliding window is set to be 10 times the foliage characteristic width, and 40 transects are employed for an image in this study [41].More detailed description of the variables in Equation ( 9) can be found in [37].Eventually the clumping index can be calculated by Equation (6-9) from digital images acquired on adjacent overcast days.

Leaf Angle Distribution and Leaf Projection Function
The leaf projection function (G) is the mean projection of unit foliage area on a plane perpendicular to the solar or view directions [34].Assuming a uniform leaf azimuth orientation distribution, the G function can be expressed as [43]: where θ is the solar or view zenith angle, θ  is the leaf inclination angle, and  =  −1 (θθ  ).
The two-parameter Beta-distribution has been evaluated to be the most accurate to describe the leaf inclination distribution function (θ  ) [44,45]: where t is 2  /, μ and ν are two parameters, and (μ, ν) is the Beta-function defined as: where Γ is the Gamma function, and  and  are calculated as: where  ̅ and σ  2 are the mean and variance of t, and σ 0 2 is the maximum variance of t, which can be calculated as: From Equation ( 10)-( 16), the leaf angle distribution (LAD) and G function can be described by two parameters,  and .Traditionally, LAD can be measured by mechanical inclinometers, while this direct method is labor-intensive and may not be feasible for tall forest canopies [46].A new developed photographic method by analyzing leveled digital camera images at different heights of the canopy allows for a rapid and accurate estimation of LAD over broadleaf canopies [46][47][48].The leveled digital images can be taken along a vertical crown profile at 2-m height increments from the crown bottom to the top [48].Leaves in the images oriented perpendicularly to the view direction of the camera are selected, and their leaf inclination angles between the leaf surface normal and the zenith direction are estimated using the angle measurement tool of an image processing software (ImageJ; http://rsbweb.nih.gov/ij/)[47].
The LAD has an impact on the regional CO2 and H2O fluxes [49], and can be affected by the trunk, stems, branches, and twigs.The spatial variability of LAD can be accounted for by dividing the canopy into multiple layers, and then the LAD of each layer and its contribution to the four components can be calculated layer by layer separately.The temporal variability of LAD can be quantified by a time-series model to extract the trend of temporal changes from LAD measurements at several times.The spatial representative of the LAD measurements may be limited in a heterogeneous forest when the tower has access to only a few species.The potential impact of light environment changes caused by the tower on the canopy structure can be reduced by the mask of canopies near the tower.In general, a minimum sample set of 75 measured leaves across the whole canopy are considered to be sufficient to characterize LAD reliably [48].The leaf inclination angle measurements are used to fit the two parameters μ and ν in Equations ( 13)-( 15), with which (θ  ) and (θ) can be derived by Equations ( 10)- (13).If no LAD field measurement data are available, then a typical LAD function for specific tree species can be assumed as the back-up approach [48].

Materials
Firstly, the different LAI extraction methods were evaluated by the 3D realistic structural scene at the YJP site.Then the proposed LAI extraction approach for sunny days was applied to extract LAI from actual images acquired at Huailai Remote Sensing Experimental Station in Beijing, China, from July to September in 2014.

Field Data Collection at the YJP Site
The YJP site in this study is part of the Boreal Ecosystem-Atmosphere Study (BOREAS) eddy covariance tower sites [14].YJP denotes young jack pine (Pinus banksiana), and is located in the BOREAS southern study area situated near Candle Lake, Saskatchewan, Canada (53.975°N, 104.650°W) as in Figure 2. The average age of the jack pine trees is 11-16 years old, and the density of the YJP site is about 4000 stems per hectare.The crown shape of the jack pine consists of a cone top and a cylinder bottom with the mean radius of 0.85 m, and the understory is composed of thin grasses, lichens, and some bearberry [50].
Three parallel transects with the equal length of 300 m, were separated by 10 m and oriented along northwest and southeast directions.The flux tower in YJP site was exactly located at the center of the middle transect.The effective LAI measurements by LAI-2000 Plant Canopy Analyzer (LI-COR, Lincoln, NE, USA) were taken along the transects every 10 m marked by the forest flags and, thus, eventually 90 readings were made within 20 min under overcast conditions.The 90° view caps were used to hide the operator from the sensor's view during the measurements by LAI-2000.The element clumping index Ω  was measured by TRAC along the transects under clear-sky conditions.The sensors were carried at about 0.1-0.2m above the forest floor with the walking speed of 1 m per three seconds.The needle-to-shoot area ratio (γ  ) was 1.43, and the mean width of foliage elements projected on a plane perpendicular to the solar direction (Ws) was 0.17 m [12].The spherical LAD (i.e., G = 0.5) was assumed, which means leaves have no preferred orientation and is often considered as a reasonable assumption for conifer shoots where no LAD measurement data are available [48].The canopy structure parameters of the YJP site are shown in Table 1.The same dataset of the YJP site has been used to evaluate the optical-based LAI measurement techniques and the improved four-scale model in previous studies [12,14,50].Table 1.The canopy structure parameters of the YJP site.Ha is the height of the lower part of the tree with no foliage; Hb is the height of the crown; R is the mean radius of the tree crowns; Ws is the mean width of foliage elements projected on a plane perpendicular to the solar direction; and G is the leaf projection function in Equation (10).

3D Reference Scene Construction
Ground truth LAI can be difficult to acquire by destructive methods in practice for tall forest canopies, and previous studies usually relied on indirect LAI measurements by optical sensors, such as LAI-2000 and TRAC, as the reference value to compare different methods [30].[12] reported that 80% accuracy can be achieved by LAI-2000 and TRAC when operated carefully, which suggests it is difficult to draw firm conclusions in method comparisons if the uncertainty of the reference values is relatively large.The 3D realistic structural scenes can account for the location, orientation, size, and shape of every individual leaf within a canopy, and thus have the advantage of knowing exactly the true values of each parameter, such as LAI, which provides the objective ground truth for the validation of various indirect LAI measurement methods [37,51].In fact, the 3D realistic structural scene based approaches have been widely used as the "surrogate truth" in the RAdiative transfer Model Intercomparison (RAMI) exercise conducted by the Joint Research Center, European Commission [52].It should be noted that the following analysis was based on synthetic images by the 3D scene instead of real images.The canopy structure parameters acquired in the field measurements at YJP site in Table 1 except γ  were used as inputs to construct the 3D scene by the 3D Max software and the Maxscript language [53].Since it was almost impossible to construct conifer forests at the needle level due to the huge number of needles, the shoot was used as the elementary unit in the 3D reference scene construction.Then the needle-to-shoot area ratio (γ  ) was set at 1 and, thus, Ω = Ω  , which was equivalent to set the forests as broadleaf trees.The constructed 3D scene has an area of 125 m 2 with 50 crowns.
For synthetic images under sunny conditions, the solar zenith angle   and solar azimuth angle   were set at −30° and 180°, respectively.Then multi-angular downward-looking synthetic images under axonometric projection in the principal plane (PP) and in the cross plane (CP) were generated at 900 × 900 pixels in 8-bit digitalization per color channel, with the view zenith angle   ranged from −50° to 50° at 10° intervals, as partly shown in Figure 3.For synthetic images under overcast conditions, the illumination was set to be isotropic and thus significant shadows cannot be seen in the synthetic images.The 3D software can determine whether each leaf was sunlit and visible at the same time by a ray tracing method.The sunlit foliage component (PT), shaded foliage component (ZT), gap fraction (GF), and true LAI can be directly extracted from the constructed 3D scene by definition as the ground truth values.The 3D reference values for clumping index (CI) at different view angles can also be derived from the 3D scene by definition in Equation ( 6), which was a more direct way than other methods such as the LX, CC or the path length distribution-based method in Section 2.3.PT was extracted from synthetic images on sunny days, GF and CI were derived from synthetic images on overcast days, while ZT was derived from all of the synthetic images because ZT was the difference between the probability of seeing the foliage and seeing the sunlit foliage.

In Situ Measurements at the Huailai Site
The Huailai Remote Sensing Experimental Station is located in Beijing, in the north of China (40.348°N, 115.783°E).The camera was fixed on the tower crane about 23 m above the ground surface, and the size of the images was 2736 × 1824 pixels.In total, 18 digital actual images were acquired under clear-sky conditions by the tower crane at nadir view (θ  = 0°) in 2014 over two different vegetation types: aspen and corn.Among them, seven images were taken over the aspen on 27 July and four, three, and four images were taken over the corn on 27 July, 18 August, and 20 September, respectively.The illumination geometry including the solar zenith angle and the solar azimuth angle for each image were determined by the location of the Huailai site and the accurate acquisition time of each image.
The heights of the aspen were measured at 35 sampled trees in the seven images by a laser distance meter, and the corn heights were measured at 20 sampled plants for each time phase.The average canopy height of the aspen was about 7.8 m on 27 July and the average corn height was about 1.6 m, 1.9 m, and 2.0 m on 27 July, 18 August and 20 September, respectively.The mean foliage diameter   for the sampled aspen leaves was 0.16 m, and   for the sample corn leaves was 0.15 m, 0.18 m, and 0.20 m.The spherical LAD was assumed for the aspen and corn as a reasonable assumption where no LAD measurement data was available [48].The size of the ground surface corresponding to each image was about 25.2 m × 16.8 m, and the four corners of the ground in the image were labeled for the LAI and the element clumping index (Ω  ) measurements.The effective LAI was measured by LAI-2000 at nine single points evenly distributed over each image after sunset of the image acquisition day.The element clumping index ( Ω  ) for the aspen was measured by TRAC along two transects under clear-sky conditions.The corn was usually considered as homogeneous canopies when the characteristic of row structure was not significant, and thus the element clumping index (Ω  ) for the corn was set as 1.The needle-to-shoot area ratio (γ  ) was set at 1 because both of the aspen and the corn were broadleaf plant species.

Scene Component Extraction
The sunlit foliage component (PT) extracted by the automated LAB2 image classification algorithm from 22 images in the principal plane and the cross plane under sunny conditions were compared with that determined by the 3D software as in Figures 4 and 5.The minimum and average overall accuracy of the LAB2 algorithm for the 22 images were 91.4% and 93.6%, respectively.The overall accuracy was the proportion of the number of pixels correctly identified as sunlit foliage component or background in the number of all pixels of the image.From visual inspection in Figures 3b and 4, the sunlit foliage were effectively extracted by the LAB2 algorithm.The slight misclassification mainly occurred on the edge of the sunlit leaf blur, which might be due to the smoothing of the penumbra effect and the light scattering and diffraction [54].Figure 5 suggests that the R 2 , root-mean-square error (RMSE) and bias for the PT extracted by the LAB2 algorithm from the 22 images were 0.98, 0.04, and −0.03, respectively.The intercept of the regression line was negative, which suggests that the PT extracted by the LAB2 algorithm was negatively biased when the value of PT was small.The slope of the regression line was larger than 1, which indicates that the LAB2 algorithm can overestimate PT when the value of PT was large.The area ratio of different scene components, including sunlit foliage (PT), shaded foliage (ZT) and gap fraction (GF) extracted by the LAB2 algorithm and the 3D reference values at different view zenith angles in PP and CP are shown in Figure 6.In PP, the sunlit foliage component (PT) reached the maximum in the hot spot direction when illumination and view directions coincide (θ  = −30°), which was due to the absence of visible shadows, but dropped quickly with the increasing of the phase angle between the directions to the sun and the camera.By contrast, a valley can be seen around the hot spot for the shaded foliage component (ZT), but the value of ZT increased significantly on both sides of the hot spot region.In CP, PT increased slightly from nadir view to large view zenith angles, while ZT increased more sharply compared with PT.In both PP and CP, the gap fraction (GF) reached the maximum at nadir view, but decreased smoothly with the increasing of the view zenith angle.Compared with PT and ZT, the gap fraction (GF) extracted by the LAB2 algorithm exhibited the best performance with an RMSE of 0.03.The shaded foliage (ZT) had the lowest accuracy with an RMSE of 0.07.

Clumping Index Estimation
The clumping index (CI) estimated by two indirect methods and the 3D reference values at different view zenith angles in PP and CP are shown in Figure 7.In both PP and CP, the CI increased with the increasing view zenith angle (  ), and reached the minimum value at nadir view, which suggests the clumping effect for such discontinuous canopies was the most significant when θ  = 0°.In PP, both the path length-based method and the LX method for CI were underestimated compared with the 3D reference values.Performances for the path length-based method and the LX method were slightly better in CP than in PP, and the accuracy of the image classification to distinguish the foliage and background in the previous step may have an impact on the CI estimation.In PP and CP, the variations of the CI by the LX method against θ  were much smoother than that by the path length-based method and the 3D reference values, which suggests that the LX method may not be very sensitive to the view zenith angle variations for the CI estimation.
The statistical results for the evaluation of the path length-based method and the LX method by the 3D reference values at different view zenith angles in PP and CP are shown in Figure 8.The clumping index (CI) was slightly underestimated by both methods with the negative bias, while the scatter points by the path length-based method were more approaching the 1:1 line compared with the LX method.The R 2 , RMSE, and Relative Error (RE) of the path length-based method were 0.86%, 0.05%, and 6.6%, while that of the LX method were 0.46%, 0.07%, and 7.8%, respectively, which suggests that the path length-based method improved the accuracy of the clumping index (CI) estimation compared with the LX method.

Comparison of LAI Retrievals with Other Methods
The retrieved LAI and the corresponding RE by four different methods at different view zenith angles in PP and CP are shown in Figure 9.The four methods include the retrieval of LAI by the sunlit foliage component under sunny conditions with the CI estimated by the path length-based method (LAI_Path) and the LX method (LAI_LX), and by the directional gap fraction model under overcast conditions (LAI_Overcast) and under sunny conditions (LAI_Sunny) with the CI estimated by the widely used LX method.In PP, the LAI_Path and the LAI_LX methods underestimated the LAI in the forward scattering directions with the view zenith angle θ  > −30°, while overestimated the LAI in the backward directions with θ  < −30°.The largest RE for the two methods in PP occured at large view zenith angles in the backward directions, which was probably due to the symmetrical characteristic of the Kuusk Hot-spot function [35].In CP, the two methods underestimated the LAI when the view zenith angle θ  ≤ 30°, while overestimated the LAI at large view zenith angles when θ  > 30° .The LAI_Overcast overestimated the LAI in PP, especially at large view zenith angles, while achieved the highest accuracy among the four methods in CP.The directional gap fraction model underestimated the LAI in both PP and CP under sunny conditions (LAI_Sunny) as in Figure 9, due to the difficulty to correctly extract the shaded foliage component from the shadowed parts of the images in the automated image classification.The misclassification of shaded foliage as background led to the overestimation of the gap fraction, and eventually resulted in the LAI underestimation by the gap fraction model.In PP, LAI_Sunny achieved the highest accuracy in the hot spot direction when the area ratio of the shaded foliage component was minimum, and the amplitude of RE increased when the view directions were away from the hot spot direction as in Figure 9a.The trend of the RE amplitude in PP was correlated with that of the shaded foliage fraction (ZT) in Figure 6a, which also increased with the increasing of the phase angle between the directions to the sun and the camera.The largest RE amplitude in PP occurred at large view zenith angle (e.g., θ  = 50°) in the forward scattering direction.In CP, the RE amplitude reached the minimum value at the nadir view, and increased with the increasing view zenith angle (θ  ) as in Figure 9b.The trend of the RE amplitude in CP was also similar to that of the shaded foliage fraction (ZT) in Figure 6b, which was symmetrical along the nadir view direction.
The evaluation results for the performances of the four different methods by 22 images in PP and CP are shown in Table 2.The LAI_Overcast achieved the highest accuracy among the four methods with an RMSE of 0.32 and an RE of 9.0%, which suggests that the directional gap fraction model was the best choice for the LAI extraction under overcast conditions.The LAI_Path method proposed in this study performed best among the other three methods using images acquired under sunny conditions, with an RMSE of 0.35 and an RE of 11.4%, which can meet the accuracy requirement (0.5; 20%) by the GCOS.The LAI_Path method performed slightly better than the LAI_LX method, which was mainly due to the more accurate estimation of the clumping index as described in Section 4.2.The LAI_Sunny method had the lowest accuracy among the four methods, with an RMSE of 1.61 and an RE of 55.9%, which was significantly outside the accuracy requirement (0.5; 20%) by the GCOS.This suggests that the directional gap fraction model was inappropriate for images acquired under clear-sky conditions due to the uncertainties in the shaded foliage extraction, and it was more accurate to extract LAI by the sunlit foliage component on sunny days.Table 2.The statistical results for the performances of four different LAI retrieval methods from 22 images in the principal plane and in the cross plane.The four methods including LAI_Path, LAI_LX, LAI_Overcast, and LAI_Sunny are the same as in Figure 9.

Applications at the Huailai Site
The actual downward-looking digital images acquired under clear-sky conditions for a corn region and for an aspen region at the Huailai site are shown in Figure 10  The comparison between retrieved LAI by 18 images acquired on sunny days with field measured true LAI at the Huailai site in 2014 is shown in Figure 11.The field measured true LAI was calculated as the ratio of the effective LAI measured by LAI-2000 and the element clumping index measured by TRAC.Two methods include the proposed LAI extraction method by the sunlit foliage component under sunny conditions with the clumping index estimated by the path length-based method (LAI_Path), and the retrieval of LAI by the directional gap fraction model under sunny conditions with the clumping index estimated by the widely used Lang and Xiang (LX) method (LAI_Sunny).It can be seen from the 1:1 line that the directional gap fraction model underestimated the LAI with the bias of −1.38 under sunny conditions (LAI_Sunny), and this was due to the misclassification of shaded foliage as background, which resulted in the overestimation of gap fraction.The proposed LAI extraction method by the sunlit foliage component (LAI_Path) was more approaching the 1:1 line with the slope closer to 1 compared with the directional gap fraction model.The proposed LAI extraction method slightly overestimated the LAI at low vegetation cover (LAI < 3), while underestimated the LAI at high vegetation cover (LAI > 5).The R 2 , RMSE, and Relative Error (RE) of the directional gap fraction model were 0.86%, 1.55%, and 36.2%, while that of the proposed LAI extraction method were 0.89%, 0.49%, and 11.6%, respectively.This suggests that the directional gap fraction model may not meet the accuracy requirement (0.5; 20%) by the GCOS on sunny days, and the proposed LAI extraction method by the sunlit foliage component improved the accuracy of the LAI estimation from (1.55; 36.2%) to (0.49; 11.6%) compared with the directional gap fraction model.

Discussion
One of the advantages for near-surface remote sensing is the high spatial resolution, which can distinguish the foliage from background and makes it possible to extract LAI continuously by near-surface imaging sensors [15,16,54].The directional gap fraction model can achieve the highest accuracy under overcast conditions (LAI_Overcast) as in Figure 9, while the performance of the traditional gap fractional model on sunny days was quite poor (LAI_Sunny), which was significantly outside the accuracy requirement (0.5; 20%) by the GCOS.The poor performance of the gap fractional model on sunny days limited the extraction of LAI by near-surface remote sensing under different illumination conditions.The proposed LAI extraction method (LAI_Path) by sunlit foliage component achieved the accuracy of (0.35; 11.4%) on sunny days, which can meet the accuracy requirement (0.5; 20%) by the GCOS.Although the accuracy of the proposed method on sunny days (LAI_Path) was slightly lower than that of (0.32; 9.0%) by the gap fraction model under overcast conditions (LAI_Overcast) partly due to the penumbra smoothing effect under sunny conditions, the proposed method improved the accuracy greatly from (1.61; 55.9%) by the gap fractional model on sunny days (LAI_Sunny).The proposed method relaxes the required illumination conditions for the extraction of LAI by near-surface remote sensing from only overcast conditions to all light conditions.The uncertainties for the non-fisheye or fisheye digital photography-based LAI extraction methods mainly come from image classification, canopy structure parameter estimation (e.g., Ω and G), and the LAI inversion model [9,47].
Illumination conditions and camera exposure settings have a crucial impact on the image classification for gap fraction estimation, and sunny illumination conditions can lead to poor gap fraction estimations due to the misclassification of the shadowed parts in the images as in previous studies [29,54].In this study, the average overestimation of the gap fraction which is, in fact, the proportion of the shaded foliage component (ZT), can be as large as 0.23 in the principal plane and 0.26 in the cross plane.The consequent underestimation of LAI is (1.61; 55.9%), which is obviously outside the accuracy requirement (0.5; 20%) by the GCOS.The underestimation of LAI on sunny days is consistent with previous studies, which strongly recommend that overcast conditions should be privileged to guarantee the image classification accuracy when feasible [9,29,30].Misclassification can occur on the edges of the leaf blur due to the light scattering and diffraction or the penumbra smoothing effect under sunny conditions, which has also been reported by [54].While the tedious and time-consuming image processing used to be considered as the main weakness for DHP by [9], the LAB2 algorithm makes it possible for the automated image processing and the accurate extraction of the sunlit foliage component with the minimum overall accuracy of 91.4%.
The Lang and Xiang (LX) method is shown to be not very sensitive to variations of the view zenith angle (  ) for the clumping index (CI) estimation in this study.The underestimation of CI by the LX method has also been reported in previous studies, which can be explained by the assumption that the foliage elements are randomly distributed within the finite length, while empty segments (no gaps) or large gaps between tree crowns will give erroneous results [37,41,55,56].The path length distribution-based method can avoid the assumption of constant path length and can effectively characterize the non-randomness within canopies, which has improved the accuracy of the CI estimation compared with the LX method.For the leaf projection function (G) and LAD estimation, previous studies have started to use multi-directional gap fraction measurements by the 5 rings of LAI-2000 or different angular sectors of the DHP image for a joint retrieval of LAI, Ω and G [9,29], while [57] argued further work was still needed to separate the effects of foliage clumping (Ω) on the estimation of G.The photographic method adopted in this study by analyzing leveled digital camera images provides an independent way for LAD and G estimation, which can avoid the influence of other canopy structure parameters, e.g., Ω and LAI [47,48].
The development of the GO model which can quantify the sunlit foliage component relaxes the required illumination conditions for the digital photography from only overcast conditions to all light conditions [33].Actually, the state of the art of forward models can significantly improve the theory and methods of ground observations.For example, the apparent clumping index proposed by [58] makes it possible for the LAI-2000 instrument to quantify the clumping effects, and the 1D bidirectional transmission model developed by [59] can overcome the traditional assumption of LAI-2000 that foliage absorbs all the radiation in the blue band, and provides a mechanism for the use of the instrument under direct sunlight conditions in the latest version of LAI-2200C.The largest overestimation of LAI by the proposed method (LAI_Path) occured at large view zenith angles on the backward side in the principal plane (e.g., θ  = −50°).This is mainly due to the underestimation of the sunlit foliage component for the current GO model in the backward directions, and this phenomenon has been reported by [53] and [60].Thus, futher development is still needed for the GO model to simulate the sunlit foliage component more accurately in the backward scattering directions.
In this study, the woody part is neglected because the shadowing effect of the trunk, stems, branches, and twigs is less of a problem for extracting sunlit foliage from downward-looking images, while the woody elements may significantly influence the extinction and gap fraction measurements for upward-pointing cameras and LAI-2000 [15,47].Additionally, the YJP site in this study does not have important understory, the canopies are relatively open and, thus, the sunlit background can be easily distinguished, while the contrast between the foliage and the background understory will have an impact on the accuracy of the image classification and the sunlit foliage extraction [38,50].It is promising to use photographs taken at different times on a sunny day with varying solar zenith angles to reduce the LAI inversion uncertainty.Finally, the digital photography method may not work due to gap saturation when canopies reach closure and the visible sunlit foliage component reaches the upper limit, especially for relatively dense tropical rainforest, which is also the limitation of other indirect methods [30].There is a little difference for the sensitivity of saturation, because the saturation is directly related to the reflectance for downward-looking sensors, such as satellite or near-surface remote sensing, while the saturation is linked to the transmittance for upward-pointing sensors such as LAI-2000.The proposed downward-looking method underestimated the LAI at high vegetation cover (LAI > 5) compared with the upward-pointing LAI-2000 in Figure 11, which suggests that downward-looking methods might be more sensitive than upward-pointing methods because the light intensity can still be decreased when penetrating the canopy even after the canopy is closed and no sky can be observed from the ground surface.

Conclusions
Recently, near-surface remote sensing using networked digital cameras provides a low-cost way to continuously monitor the vegetation dynamics over the regional to continental scale.However, current indirect LAI measurements by the directional gap fraction model can only use images acquired under diffuse sky conditions, while tower-based webcam images can also be acquired under clear-sky conditions with the direct sunlight.The main challenge for the images on sunny days is that it is difficult to discriminate the shaded foliage and shaded background from the shadows.A new LAI extraction method by the sunlit foliage component from downward-looking digital photography under clear-sky conditions is proposed in this study.This method extracts the sunlit foliage component by the automated LAB2 image classification algorithm, estimates the clumping index by a path length distribution-based method, quantifies the LAD and G function by leveled digital images, and eventually obtains the LAI by introducing a GO model which can quantify the sunlit foliage fraction.The proposed method was evaluated at the YJP site, Canada, by the 3D realistic structural scene constructed based on the field measurements.Then the proposed LAI extraction approach for sunny days was applied to extract LAI from actual images acquired at Huailai Remote Sensing Experimental

Figure 1 .
Figure 1.The flow chart of procedures to extract LAI from downward-looking digital photography under clear-sky conditions.

Figure 2 .
Figure 2. The young jack pine (YJP) site located in the BOREAS southern study area situated near Candle Lake, Saskatchewan, Canada (53.975°N, 104.650°W).The projection is UTM 13 North, WGS84.The background image is a true-color composite image from Landsat8/TM that was acquired on 8 August 2014.

Figure 3 .
Figure 3.The downward-looking synthetic images from constructed 3D scene of the YJP site under clear-sky conditions.The components in bright green, bright brown, dark green, and dark brown correspond to sunlit foliage, sunlit background, shaded foliage and shaded background, respectively.(a-c) correspond to synthetic images generated in the principal plane with the view zenith angle θ  of 30°, 0° and −30°, respectively.The solar zenith angle θ  and solar azimuth angle φ  are −30° and 180°, respectively.

Figure 4 .
Figure 4.The classification results of sunlit foliage (in black) and background (in white) by the automated LAB2 algorithm (a) and by the 3D software (b).The original RGB image was acquired at nadir view (θ  = 0°) as in Figure 3b.

Figure 5 .
Figure 5.Comparison between the area ratio of the sunlit foliage component (PT) extracted by the LAB2 algorithm (on the vertical axis) and the 3D reference values (on the horizontal axis) from 22 images in the principal plane and in the cross plane.

Figure 6 .
Figure 6.The area ratio of the sunlit foliage component (PT), shaded foliage component (ZT) and gap fraction (GF) extracted by the LAB2 algorithm and the 3D reference values in the principal plane (a) and in the cross plane (b) with the solar zenith angle (SZA) of −30°.

Figure 7 .
Figure 7.The clumping index (CI) estimated by the path length distribution-based method in this study (CI_Path) and the LX method (CI_LX), compared with the 3D reference values (CI_3D) in the principal plane (a) and in the cross plane (b).

Figure 8 .
Figure 8.Comparison between the clumping index (CI) estimated by the path length distribution-based method (CI_Path) or the LX method (CI_LX) with the 3D reference values (CI_3D) from 22 images in the principal plane and in the cross plane.

Figure 9 .
Figure 9.The retrieved LAI and the corresponding Relative Error (RE) by four different methods in the principal plane (a) and in the cross plane (b).The LAI_True is the true LAI set in the 3D realistic structural scene, which serves as the ground reference value; LAI_Path and LAI_LX are the retrieval of LAI by the sunlit foliage component under sunny conditions with the CI estimated by the path length-based method and the LX method, respectively; LAI_Overcast and LAI_Sunny are the retrieval of LAI by the directional gap fraction model under overcast conditions and under sunny conditions, respectively.

Figure 10 .
Figure 10.The actual downward-looking digital images acquired under clear-sky conditions for a corn region on 20 September 2014 (a) and for an aspen region on 27 July 2014 (b) at the Huailai site.Both of the two images were acquired at nadir view ( θ  = 0°).The classification results of sunlit foliage (in black) and background (in white) by the automated LAB2 algorithm for the corn region and the aspen region are shown in (c,d), respectively.

Figure 11 .
Figure 11.Comparison between retrieved LAI and the corresponding Relative Error (RE) by 18 images acquired on sunny days with field measured true LAI by LAI-2000 and TRAC at the Huailai site in 2014.LAI_Path is the retrieval of LAI by the sunlit foliage component under sunny conditions with the clumping index estimated by the path length-based method, and LAI_Sunny is the retrieval of LAI by the directional gap fraction model under sunny conditions with the clumping index estimated by the widely used Lang and Xiang (LX) method.