Crown level clumping in Norway spruce from terrestrial laser scanning measurements

The clumping of coniferous needles into shoots is widely acknowledged as a structural feature that cannot be ignored in radiation regime models of coniferous forests. However, higher level clumping, i.e. the aggregation of leaves and shoots in tree crowns and forest stands, is still rarely accounted for in the models. Clumping reduces the light interception of and increases the light penetration depth in a plant stand. To improve forest radiation regime models with respect to this forest structural parameter, we propose a method that can quantify clumping at different hierarchical levels by estimating the silhouette to total area ratio from point clouds acquired by laser scanners. Our method is based on estimating attenuation coefficients in a voxel grid, and subsequently computing the total leaf area and spherically averaged silhouette area of a tree crown or forest stand. We tested our method with empirical data in young Norway spruce trees, where we compared leaf area and silhouette area to destructive and photogrammetric reference measurements. The accuracy of leaf area estimates depended strongly on the voxel size, with voxel sizes below 10 cm side length exhibiting up to 100% higher estimates than the reference leaf area, and large voxels with 90 cm side length being closest to the reference measurements due to crown clumping. The silhouette area estimates varied less with voxel size and were slightly higher than the reference estimates. We analyzed possible error sources and point out ways to improve the measurements of leaf and silhouette area for conifer trees using laser scanning data.


Introduction
The structure of a forest canopy influences the absorption and scattering of solar irradiation, and thus plays a key role in determining primary production, and hence in forest growth and carbon sequestration. Forest canopy structure is typically described through the density, orientation and positions of leaves within the canopy. Typically, these three aspects of forest canopy structure are parameterized by the total leaf area of a stand, the leaf angle distribution, and the spatial dispersion patterns of leaves. The total leaf area is often expressed as a relative measure, either as leaf area index at the stand level, or as a leaf area density in reference volumes. The spatial dispersion pattern is mostly called clumping due to the typically aggregated patterns that have been found in forest canopies. Theories as well as measurement methods have been built on these three structural characteristics of foliage, usually linking foliage structure to the radiation regime. Early measurement methods focused principally on agricultural crops, of which one of the best known is the point quadrat method (Wilson, 1960), which estimates the leaf area index and the radiation extinction coefficient in plant stands. The extinction coefficient originates from the Beer-Lambert law, which was adapted to plant stands by adding a term that considers the orientation of scattering elements, which is more suitable to describe leaves (Ross, 1981). Nilson (1971) built the foundation of many following measurement techniques by developing a theoretical framework which enables the retrieval of leaf area index from measurements of angular gaps in canopies. Gap-based methods have since become the standard for indirect measurements of leaf area index in forest stands. A variety of techniques is available, ranging from hemispheric photography to specialized instruments such as the LAI-2200 Plant Canopy Analyzer (Fournier and Hall, 2017).
In the past decades, advances in Light Detection And Ranging (LiDAR) technology and methodology have enabled us to collect threedimensional data at an unprecedented rate and resolution. The different LiDAR scanning platforms offer different relative advantages, which, besides different viewing geometry are associated with occlusion, or the penetration depth of LiDAR beams into canopies (Morsdorf et al., 2018).
In a nutshell, choosing LiDAR platforms faces a trade-off between high resolution (best achieved with terrestrial platforms) and large area coverage (airborne platforms). Mobile and UAV platforms offer a compromise between resolution and large area coverage. The resulting point cloud data contains detailed geometric information about forest canopies, with the main disadvantage being the lack of topological information. That is, a point cloud is an unordered list of coordinates that does not contain information on the objects present in a scene. Nevertheless, a plethora of methods for processing forest point clouds has been established that either aim at reconstructing topology information (e.g. Atkins et al., 2018;Hackenberg et al., 2015;Raumonen et al., 2013) or directly use point clouds to estimate forest structural characteristics. All three foliage characteristics, i.e. leaf area, leaf angle distribution, and clumping, have been studied with various methods. Leaf area has been estimated using either gap-based methods (e.g. Zhao et al., 2015) or voxel-based methods (e.g. Béland et al., 2014;Pimont et al., 2018). Leaf angle distribution can be estimated simultaneously with leaf area for entire stands by simultaneous retrieval of the orientation parameter G and the leaf area of tree crowns (Zhao et al., 2015), or, if leaves are resolved in the point cloud, by manual angle measurements of individual leaves from the point cloud (Béland et al., 2011;Zhu et al., 2018).
Clumping has recently received increased attention in the terrestrial laser scanning (TLS) modeling community. Currently, two approaches have been applied to quantify clumping from laser scanning data, the gap size distribution (Li et al., 2017;Zhu et al., 2018) and the path length distribution (Chen et al., 2018). Both approaches share the basic principle of analyzing gap size realizations and compare them to the gap size distribution of a Poisson canopy. The aim of clumping quantification in this sense has been to produce a clumping index, i.e. a correction factor used to convert the true leaf area index to an effective leaf area index with random foliage dispersion that is used in the modified Beer-Lambert law.
However, clumping is more than a correction factor, as it is an important forest structural characteristic which influences the fraction of absorbed and scattered shortwave radiation of forests, the penetration depth of solar irradiation into a stand and thus the light availability to both the lower parts of tree crowns and the understory vegetation (Kim et al., 2011). It has become commonplace to correct for shoot clumping in conifer forests. Shoot clumping was originally quantified by Oker--Blom and Smolander (1988) using the spherically averaged silhouette to total area ratio (STAR). Subsequently, shoots were considered as the basic scattering elements in coniferous forests, and gap-or transmission-based estimates of leaf area were corrected by the clumping index 4STAR. Forest radiation regime models, based on e.g., the spectral invariants theory (Knyazikhin et al., 1998), can distinguish between wavelength-independent (i.e., forest structural) parameters and wavelength-dependent (i.e., spectral) parameters. The influence of forest structure on its radiation regime can be described by a single parameter, the photon recollision probability (Smolander and Stenberg, 2005). The concept of photon recollision probability has been linked to shoot STAR, indicating the hierarchical nature of scattering in conifer canopies (Smolander and Stenberg, 2003). The photon recollision probability was also used to create the PARAS reflectance model (Rautiainen and Stenberg, 2005). Inversion studies of PARAS, however, showed a tendency to overestimate leaf area index (Heiskanen et al., 2011;Varvia et al., 2018). Both studies used inversion models that did not consider clumping, which likely caused PARAS to overestimate reflectance, and therefore, in the inversion, lead to an overestimate of leaf area index. Schraik et al. (2019) showed that particularly for forest stands with low leaf area index, the clumping index can reach levels lower than what has been observed for shoot clumping alone. Simulations on Scots pine trees indicate that there may be a strong effect of foliage clumping at hierarchical levels above shoot scale (Stenberg et al., 2014). However, this effect of higher level clumping has not yet been empirically validated. If clumping at higher levels than the shoot can be shown to be non-random, its incorporation in forest radiation regime models will likely improve the models' ability to accurately describe the relationship between forest structure, biochemical properties, and optical properties of vegetation elements and forest stands. In addition, knowledge on higher level clumping may help to improve estimates of leaf area by relaxing the assumption of random distribution of leaves in a canopy.
In this article, we set out to develop a method for estimating clumping using LiDAR point cloud data. We study the effects of clumping within individual tree crowns using empirical TLS point clouds. The method we propose is based on estimating the two components of STAR, i.e., the total leaf area and the spherically averaged silhouette area. To validate our method, we compare our TLS-based method for estimating tree crown silhouette and total leaf area to photographic and destructive measurements. Secondly, we analyze how technical aspects, such as voxel size and scan resolution, influence the prediction accuracy of silhouette and total areas. Finally, we use our method to demonstrate how STAR of individual tree crowns varies in young Norway spruces.

STAR as a clumping descriptor
According to one of Cauchy's theorems, the spherically averaged silhouette to total (surface) area ratio of a convex body is 1/4 (Lang, 1991). STAR of a shoot is typically smaller than 1/4, indicating mutual shading of needles within a shoot. Shoot STAR is the ratio of the shoot's average silhouette area SA to the total surface area TA of all needles within the shoot (Oker-Blom and Smolander, 1988). Since the needles are assumed convex, their spherically averaged projection area is TA/4. Consequently, 4STAR n = SA/(TA /4) is the ratio of shoot to needle silhouette area, thus quantifying needle clumping in shoots (Nilson et al., 1999;Smolander and Stenberg, 2003;Stenberg et al., 1994).
STAR can be applied to higher hierarchical levels using porous envelopes to describe the individual elements (Silva et al., 2008). Stenberg et al. (2014) used this multiscale approach to decompose the total clumping index Γ into several hierarchical levels: clumping of needles in shoots Γ n , of shoots in crowns Γ s , and of crowns in stands Γ stand . The product of these clumping indices yields the total clumping index Γ = Γ stand Γ s Γ n .
(1) Stenberg et al. (2014) showed that the clumping of needles into crowns can be formulated as Γ crown = Γ s Γ n = 4STAR crown , where STAR crown is defined as spherically averaged crown silhouette area to total needle area in a crown. In this study, we will use STAR crown to quantify crown-level clumping.

Overview of the method
In previous studies, measurements of STAR shoot were aimed at measuring its two characteristics, the silhouette area and the total needle area, using photogrammetric and destructive methods, respectively. A measurement of STAR crown can be made in a similar fashion, i.e. measuring or estimating the crown silhouette area and the total needle area of a crown. This approach is independent of data collection frameworks, which can be done from three-dimensional models or from TLS point clouds. The former framework appears more attractive due to the topology information associated with the three-dimensional structure of tree crowns. However, three-dimensional reconstruction methods have mostly focused on the woody structures of trees, while reconstruction of leaves remains difficult (e.g. Åkerblom et al., 2018;Hackenberg et al., 2015;Raumonen et al., 2013). Therefore, we see the applicability of three-dimensional models currently limited to simulated data.
TLS point clouds, on the other hand, provide no topology information, and are thus not able to distinguish between different hierarchical levels of clumping without preprocessing the data to segment individual elements such as tree crowns in stand level data. Despite the lack of topology information, TLS point clouds appear to be, at the time of writing this article, the most suitable basis for estimating the silhouette and total area of tree crowns.
The method we propose is based on the spatially explicit estimation of attenuation coefficients in a voxel grid that covers the entire point cloud. Using the attenuation coefficient, it is possible to estimate the leaf area density as well as the projection area within the same theoretical framework. This physical modeling framework enables comprehensive analysis of the point cloud data with respect to the measurement principles of TLS, and allows the assumption that modeling errors are shared between the two components of STAR as far as their common foundation is concerned, i.e. the attenuation coefficient. Pimont et al. (2018) showed that the modified contact frequency (MCF) is the best predictor for the attenuation coefficient, given the assumptions of random leaf distribution inside the voxels and an infinitesimal TLS beam footprint. The MCF is based on the point quadrats method (Wilson, 1960), which assumes a plant stand is probed an infinite number of times and the number of leaf contacts for each probe is recorded. Based on this contact frequency, one gets the attenuation coefficient λ as

Estimating total leaf area
with the number of contacts C for probe j over the number of probings N with probe length δ. To apply this method to TLS point clouds, the differences in measurement principles have to be accounted for. The main differences are that in TLS, the "probe" decays after the first contact (for single return instruments), the length of the probe is variable, the vegetation elements have a finite size, and the number of probings can be low for some voxels due to occlusion or small voxel size, for example. Béland et al. (2011) introduced the MCF to account for the decay after the first contact by including the length of the path that was actually explored by the TLS beam, i.e. the free path z, as which models the contact frequency as an indicator function that has the value one if the beam is intercepted in a voxel, i.e. the free path z is shorter than the theoretical path length δ of the ray through the voxel. This modification accounts for unequal path lengths and the decay of the probe after the first contact. Béland et al. (2014) modified the estimator of λ to consider the finite size of vegetation elements with respect to the voxel cross section size. The MCF was further corrected in Pimont et al. (2018) for theoretical biases due to a finite number of beams that enter each voxel. Based on these corrections, we calculated the attenuation coefficient per voxel as with the ratio of intercepted to transmitted beams I, the mean element size-corrected free path length z e , and the element-size corrected fraction of intercepted beams 1 zj<δj . The first term of the right handside corresponds to the modified contact frequency, and the second term is a correction term that becomes negligible with increased exploration of the voxel volume. The approach we used to obtain attenuation coefficient is described in more detail in Fig. 12 of Pimont et al. (2018). The leaf area density for each voxel can be estimated if the average projection area of unit foliage area is known, which is defined in the Gfunction (Ross, 1981). The G-function quantifies the average projection area of a unit leaf area in a zenith angle θ. In TLS-based studies of leaf area density, the value of G is often assumed constant (independent of zenith angle) at G = 0.5, corresponding to its spherical average in flat leaves (Miller, 1967). The effective one-sided (i.e. half of total) leaf area density is calculated as LAD e = λ/G. Provided that the voxel size is chosen so that the spatial distribution of leaf elements within voxel is random, LAD e equals LAD. For conifers, however, the clumping of needles in a shoot poses a challenge, because it is practically impossible to make the voxel size small enough to ensure random distribution of needles within a voxel. Additionally, due to the non-zero footprint size, the smallest unit that is resolved in the point cloud is the individual shoot (Ma et al., 2018). One can therefore consider LAD e as an effective shoot silhouette area density, which corresponds to the hemisurface needle area LAD as with the factor β quantifying the effect of the effective beam footprint , and 4STAR shoot accounting for clumping of needles into shoots. Our purpose in the current pilot study is to compare LAD e to reference measurements of LAD, i.e. we do not apply correction factor α to our results. This is in order to maintain comparability with existing studies. In practice, LAD could be estimated if STAR shoot and β are known. STAR shoot can be estimated with photogrammetric and destructive methods, and values for different species have been published in literature (Thérézien et al., 2007). Determination of β and its dependence on scan settings is a subject of a separate follow-up study. Finally, the one-sided (half of total) leaf area (LA) can be obtained by summing the leaf area density over the voxel grid.

Estimating silhouette area
The directional projection area of a voxel grid can be estimated by casting rays through the voxel grid and calculating the absorbed energy using the Beer-Lambert law. The parallel rays are equally spaced and form a synthetic image with a defined pixel size, and the pixel value corresponding to the fraction of absorbed energy along the ray, or a shadow fraction of the pixel area. The spherically averaged silhouette area SA is estimated as the projection area P of the voxel grid, which depends on λ, averaged over the spherical directions (ϕ, θ) as Finally, the silhouette to total area ratio at crown level is calculated as STAR crown = SA/(2LA).

Materials and methods
We collected 15 young Norway spruce (Picea abies L.) trees with height below 4 m to test our method for estimating STAR crown and to compare the estimates for total leaf area and directional projection area to results from destructive leaf area measurements and photogrammetric measurements, respectively. The trees were taken from a managed forest area in Southern Finland (near 60 • 23 ′ N, 24 • 37 ′ E). We selected 5 trees from a plantation area (tree numbers 1-5), and 10 trees from the understory of mature stands (tree numbers 6-15). The trees are illustrated in Fig. 1, along with data on tree height and branch (fresh) biomass.
All trees were obtained from a small forest area (within about 150 m from each other), hence the growing conditions, except for light availability, were equal for all trees. The trees were manually harvested and brought to an outdoor measurement facility at Aalto University's Otaniemi campus. We tied water-soaked paper towels to the cutting area of the trees to avoid drying of plant tissue during transport and storage. The trees first underwent outdoor measurements, where we collected TLS data and photographed the trees. Subsequently, the trees were clipped into branches, and measured in the laboratory to determine leaf area and biomass. We ensured that the branches were scanned and clipped within 30 hours of harvesting, and then stored in a refrigerator for less than 48 hours.

Measurement layout for terrestrial laser scanning and photographs
We used the TLS device Leica P40 Scan Station to scan the trees from six directions, and took simultaneous photographs. The laser has a wavelength of 1550 nm with a beam divergence of 0.23 mrad and a beam size of 3.5 mm upon exiting the instrument. The P40 Scan Station records a single-return point cloud based on waveform digitization.
Each sample tree was fixed between two metal bars that were driven into the ground to keep it stable for the duration of the measurements. We scanned each tree from six directions at 10 m distance, with a spacing of 60 • between scan locations (Fig. 2). We scanned at two resolutions, 3.1 mm at 10 m (0.018 • ) and 1.6 mm at 10 m (0.009 • ). For coregistration of the lidar scans, five circular black and white targets (Leica 4.5′′ B&W targets) were placed evenly throughout all directions (centered around the tree mount) at about 15 m distance. The coregistration of lidar scans, and a manual deletion of points on the ground or the metal mount were carried out in Leica Cyclone (version 9.4.0). Besides the manual cleaning of the tree mount and ground points we did not manipulate the point cloud any further, i.e. no filters were applied.
The TLS point clouds, as exported by Leica Cyclone in the e57 format, contain information on return pulses, and the origin of the points. Since Cyclone version 9.4, pulses which do not trigger a return (hereafter called empty pulses) are no longer included in exported point clouds. Empty pulses contain important information that is needed for estimating voxel-based attenuation (e.g. Soma et al., 2018). We recovered the empty pulses by converting the point cloud into spherical coordinates and inserting empty pulses with dummy range value (1000 m) into the point cloud (see Appendix A for details).
Photographs were taken concurrently with the scans by holding the camera with the lens positioned just above the scanner's oscillating mirror, thus reaching a quasi-coincident viewing geometry. Before the photographs were taken we set up a white background frame (a wooden frame with attached white fabric spanning approximately 2.4 x 2.8 m) to enable separation of the tree from its background. We used a Sony A7R camera with a 28 mm lens.

Destructive measurements of leaf area
We measured 20 shoots per tree, which we randomly picked from the bulk of branches that consisted of all shoots with twig diameter less than 2 cm. This was done by setting a random timer that went off between 5 and 30 seconds while selecting one shoot after another in rapid succession. The shoot that was selected when the timer went off was added to the sample. This process was repeated for all shoots. We weighed each shoot, plucked its needles and scanned them on a flatbed scanner, then weighed the needles separately. The twigs were weighed and scanned in the same manner. The needle shape was approximated by a parallelepiped (Homolová et al., 2012). Therefore we measured the major and minor axis diameter of two needles per shoot with a caliper. The flatbed scan images were processed with Python and Scikit-Image (van der Walt et al., 2014) to extract the number and average length of needles.
The resulting data comprised for each shoot its total mass, needle and twig mass, needle major and minor axis diameters, the average needle length, the number of needles, and the projected area of needles in a shoot. In addition to these shoot level data, the mass of all shoots with twig diameter less than 2 cm was recorded at tree level. We calculated the average shoot mass to total leaf surface area ratio and upscaled to the tree level using the mass of all shoots, thus arriving at the total (twosided) leaf surface area per tree. Hereafter, we will use the one-sided (i.e. half of total) leaf surface area for easier comparison with results published in other articles, except for the calculation of STAR, which by definition requires the two-sided (i.e. total) leaf area.

Photogrammetric projection area estimates
The photographs of trees were processed to reference estimates of projection area in six horizontal directions per tree. We manually cropped the photographs at the edges of the background frame, allowing some extra space at the lateral and top edges, but not at the bottom edge. The RGB images were converted to greyscale and thresholded using Otsu's method (Otsu, 1979). In about 10% of the images, where photos were taken in direction of the sun, there were shadows from the background frame that were thresholded along with the tree. In these cases, we reduced the threshold value from Otsu's method by 20% to avoid these shadows. To remove any other object from the image, we masked out any thresholded regions that were connected to the lateral and top edges. We determined the projected pixel size using the camera's sensor size and resolution, and the focal length from the specifications of the camera and lens, respectively. In the calculation of the projection area, we assumed the photographs to be taken spatially coinicident with the TLS scans, i.e. the distance between the image and the midpoint of the tree was about 10 m. We further assumed the lens distortion was negligible because the trees were centered in the image and were within a narrow fraction of the field of view of the camera.

Estimates for leaf and projection area from TLS data and comparison to reference estimates
We estimated the attenuation coefficient in a voxel grid for each tree using the modified contact frequency as presented in Pimont et al. (2018, Fig . 12 therein) and reviewed in Section 2.2.2 of our study (Eq. 4), using the ray tracing algorithm by Amanatides and Woo (1987). From the attenuation coefficient we calculated the leaf area density per voxel as and then arrived at the one-sided (half of total) leaf area per tree by summing over the voxel grid The value G = 0.574 is derived from Stenberg (2006), as in our point cloud data the shoot is the smallest resolved element. This higher G value accounts for the fact that a shoot is better resembled by a cylinder than a flat object, and hence the average projection area per unit leaf area is slightly increased. The projection area was estimated by generating synthetic images as explained in Section 2.2.3, in six equally spaced horizontal directions that corresponded to the directions from which the TLS scans and photographs were taken. The resolution of the synthetic images was 2 cm, less than half of the smallest voxel size, to avoid potential aliasing effects while maintaining computational efficiency. For comparison with the photogrammetric data, we calculated the mean of the projection area for each tree from both TLS-based and photogrammetric measurements.
For each tree, we estimated leaf area and average projection area for voxel sizes between 5 and 40 cm in increments of 5 cm, and for 50 to 90 cm in increments of 10 cm, i.e. 12 different voxel sizes were tested, corresponding to the voxel sizes used in Soma et al. (2018). We scanned each tree at two scan resolutions, at 3.1 mm at 10 m, and at 1.6 mm at 10 m. Point clouds of each resolution (0.018 • and 0.009 • ) data were used separately, and in addition we downsampled each point cloud to half resolution. Thus, 4x12 sets of results were generated (one for each scan resolution and one for each voxel size).
For each voxel size and scan resolution, we compared the leaf area from our proposed method with the destructive measurement data, and the horizontally averaged projection area with the photogrammetric data. In the comparison of photogrammetric data, we had to exclude tree number 13 because its size was too large for the background frame, hence it was not possible to estimate the tree's projection area from the photographs. To quantify differences between point cloud estimates and reference measurements, we used the relative mean deviation or bias and the relative root mean square error Table 1 Relative difference between TLS and reference estimates for total leaf area, in percent, for all voxel sizes and scan resolutions. The shorthand DS in the resolutions refers to "downsampled".
with the reference measurements y i for tree number i, and their mean y, the estimates from point cloud data ŷ i , and the number of trees n.

Estimates of STAR crown
Finally, STAR crown was computed for each tree, using the leaf area and spherically averaged silhouette area estimates from TLS data. For the estimation of spherically averaged silhouette area, we computed synthetic images as explained above, but now for 72 directions that corresponded to nodes of product Gaussian quadrature (Atkinson, 1982). This allowed computation of spherically averaged silhouette area using numerical integration. Total leaf area was computed from the estimated one-sided leaf area by multiplying with 2. Similarly as leaf and silhouette areas, also STAR was estimated separately for 12 different voxels sizes and 4 scan resolutions.

Leaf area estimation
Overall, the leaf area estimates produced by the modified contact frequency were higher than the destructive reference estimates (Table 1, Fig. 3). The degree of overestimation by MCF depended strongly on the voxel size, and only marginally on the scan resolution except for the smallest (5 cm) voxel size, for which using the lowest scan resolution deteriorated the accuracy (Table 1) due to insufficient exploration of the voxel volume. The effect of scan resolution has been described in Pimont et al. (2018), who showed that a low beam number per voxel can exhibit a strong positive error on the attenuation coefficient. For the highest resolution, we observed an overestimation of leaf area between 86% in the smallest, and 1% in the largest voxels (Fig. 3).
The deviations between MCF and reference leaf area estimates are likely due to convolution of several error sources, which cannot be separated from each other using the material and methods used in this study. The error sources are caused mainly by the finite footprint size of the TLS instrument, and within-voxel clumping. We further distinguish within-voxel clumping into clumping of needles in shoots, and clumping of shoots in a voxel i.e. nonrandom distribution of shoots inside a voxel. The finite size of the laser footprint causes an overestimation of attenuation, as a pulse needs not be centered on a leaf to produce a return. The gaps within conifer shoots are often smaller than the beam footprint, causing the within-shoot gaps to be underestimated for parts of the shoot above the return threshold, and therefore the shoot silhouette area, and thus the attenuation coefficient, is overestimated. Visual examination of point clouds supported this hypothesis, because the Norway spruce shoots were resolved as almost fully opaque cylinder-like shapes in the point clouds. On the other hand, using attenuation coefficient for conifers without correcting for shoot clumping leads to an underestimate of leaf area. Therefore, there are two errors with mutually opposite effects appearing at the scale of the individual vegetation element, and thus are present for all voxel sizes.
With increasing voxel size we observed a decreasing bias of leaf area, which can be attributed mainly to clumping of shoots within voxels by incorporation of empty voxels between branches. In addition, larger voxel size, depending on the origin of the voxel grid, increases clumping in voxels at the outer crown boundary by including empty space outside the crown. Therefore, we argue that the low bias for the largest voxel size did not produce the best estimation accuracy but rather indicated that the clumping of shoots in voxels compensated for the overestimation caused by finite footprint size.
There are other possible error sources which may have played a role in our analysis. Namely, we did not separate wood and leaf returns in the point clouds, and we did not explicitly account for the multi-view geometry in our measurement setup (Pimont et al., 2019). A wood-leaf separation may have slightly improved our results, because our reference estimates did not include wood areas. However, the trees in our sample were small and woody elements were only clearly visible at the base of the stem. The trees were mounted so that the crown base was almost immediately above the tree stand, which was manually removed from the point clouds, and therefore only a very small portion of the stem was visible. Accounting for the multi-view geometry can be done by assigning a direction-dependent value for the G-function (Pimont et al., 2019). In our data, all scans varied only in azimuthal direction, for which the G-function is commonly assumed invariant, and also estimating a leaf angle distribution from conifer shoots is uncertain due to lacking resolution of needles by the TLS instrument. The G-function, besides sub-voxel clumping and beam size effects, is one of the most important factors in measuring leaf area density of trees, and subsequently STAR. However, the value for G cannot be directly estimated from our data, because it is multiplied with biases from sub-voxel clumping and beam size, among others. Resolving these biases is the focus of future work.
The tendency of the modified contact frequency (in the form proposed by Pimont et al. (2018)) to overestimate attenuation and thus leaf area has already been demonstrated in Soma et al. (2018). Soma et al. (2018) aimed at introducing empirical correction factors to account for effective beam size. The correction factors were calculated from the ratio of destructive to TLS-based leaf area α = LA ref /LA TLS . They reported values for α ranging from about 0.5 at 9 cm voxel size, to about 2.0 at 70 cm voxels. We noticed a similar relationship between the ratio α and voxel size, but in our data α ranged from 0.56 for 5 cm to 1.05 for 90 cm voxel size. Our results for the range of α are likely due to species differences, and differences in instruments and measurement setups.
Finally, we will discuss the two extremes, tree 5 and tree 6, which appeared as outliers in Fig. 3B. These extreme cases highlight practical examples of the above discussed errors in estimating leaf area through a voxel-based approach, and demonstrate how a potential violation of the assumption that LiDAR beams have an infinitesimal footprint can lead to both under-and overestimation of leaf area density. For most trees, the 90 cm voxel size came closest to the reference estimates, and only in one tree it was the smallest voxel size (5 cm) which was closest to the destructive leaf area estimate. The latter tree was tree 5, which was among the smallest trees in our sample, but at the same time it had a relatively high biomass of 1.42 kg compared to the similarly sized trees 2 and 9, which had 0.43 and 0.48 kg, respectively (Fig. 1).
Tree 5 had one of the densest crowns, which likely caused occlusion of TLS beams towards the inner parts of the crown. Due to this limited exploration of the inner crown volume, the attenuation of the entire tree (i.e. for voxel sizes where the entire crown is contained in a single voxel) may appear similar to trees 2 and 9, despite the much higher density. For smaller voxel sizes, the inner parts of the crown were not explored, which likely led to an underestimation of insufficiently explored voxels. Such an error can be partly compensated by classifying insufficiently explored voxels and using neighboring voxels to estimate their leaf area density, however, unexplored voxels were very few in our study and we did not investigate such a correction further.
On the other extreme is tree 6, which was most overestimated by the TLS-based leaf area. Tree 6 had a biomass of about 0.7 kg, but occupied a relatively large volume (Fig. 1). With this relatively low biomass spread over a large volume, the leaf area density per voxel was still high enough to produce LiDAR returns at a similar rate as a very dense tree. In this case, the leaf area was more evenly spread out over the voxel volumes, and appeared to be indistinguishable from an actually dense tree because the leaf area density was high enough to trigger TLS returns, thus exhibiting an extreme example of error caused by the violation of the assumption of infinitesimal beam size in the MCF.

Silhouette area estimation
Similarly to the leaf area, the silhouette area predicted by the TLSbased method was higher than the photogrammetric reference measurements. Compared to the leaf area estimates, the magnitude of overestimation of silhouette area was smaller and did not depend as strongly on the voxel size. The relative bias ranged from -22% at 5 cm voxel size, to -28% 90 cm voxel size (Fig. 4, Table 2). Similar to the leaf area, the smallest voxel size had an increased bias at resolutions of 0.018 • and higher.
We expected that our proposed method would overestimate silhouette area to some extent, as the underlying attenuation coefficient was also overestimated. In addition, silhouette area is generally overestimated by TLS methods if no correction for the effective beam size is considered. Such a correction can be either done on the attenuation coefficient, as demonstrated in Soma et al. (2018), or by using the intensity of each pulse to empirically estimate the relative gap fraction of each beam (Hancock et al., 2014). For practical purposes, the former approach would require more empirical study to produce reliable correction factors that take into account potential differences between tree species, measurement setups, and instrument specifications. The latter approach has not yet been adapted to the framework of the modified contact frequency, but would have the advantage of not

Table 2
Relative difference between TLS and reference estimates for silhouette area (SA), in percent, for all voxel sizes and scan resolutions. The shorthand DS in the resolutions refers to "downsampled".

Silhouette area
High ( requiring any calibration to species but only to instrument specifications. Besides the overall bias of the attenuation coefficient, we found that the dependence of silhouette area estimates on voxel size was not as straightforward as with leaf area index. In our results, there was a generally increasing trend of silhouette area with increasing voxel size. The largest deviation between trees appeared at the smallest voxel size, where for larger trees the silhouette area showed a decreasing trend between 5 and 15 cm voxel size, followed by a small but increasing trend (Fig. 5). The smallest trees, on the other hand, had an increasing trend throughout the voxel size range, except for a premature peak around 30 to 50 cm voxel size. This premature peak for small trees was likely caused by a processing artifact due to the placement of the voxel grid relative to the crown position. At voxel sizes of 30-50 cm there was a large fraction of voxels that were located at the crown boundaries, and thus overlapping only partly with the crown, which caused overestimation of silhouette area. At larger voxel sizes the fraction of these "boundary voxels" was lower, due to crown dimensions matching the multiple of voxel side length (Fig. 5).
The potential sources of bias in estimation of the attenuation coefficient have already been discussed in Section 4.1. In addition to the biases present in the attenuation coefficient, the silhouette areas were estimated through a second ray tracing step which used the Beer-Lambert law to estimate the absorbed energy (relative to unity) along the paths through the voxel grid. This second ray tracing step propagated, and possibly amplified errors present in the attenuation coefficients. An amplification of errors could have occured particularly in voxels located at the crown boundaries. In such boundary voxels, all the  The median of gradients of the means is the angular resolution δζ; The mean of index 0 is the starting point ζ 0 ; end Construct a grid of the size of the maximum index values from φ and α angles; Based on the grid index (u, v), calculate the nominal direction as φ nom = φ 0 + uδφ (α nom is calculated analogously with α and v); Assign measured angles to the grid based on their indices; Rearrange the grid to get a collection of vectors [φ, α] T ; Append range values of measured points; For empty points, append a dummy range value (e.g. 1000 m); Transform the vector [φ, α, r] back to Cartesian coordinates; Transform to the global point cloud coordinate system; Assign the origin of each point, i.e. the scanner position; end Algorithm 1. Recovery of empty pulses from a return-only point cloud crown volume is concentrated on one side of the voxel, therefore being highly clumped. By an overestimated attenuation coefficient for such boundary voxels, the empty space outside the crown but within the voxel contributes to the silhouette area, thus possibly causing the slightly increasing trend in silhouette area with increasing voxel size. However, the variation in silhouette area between the extreme voxel sizes was small (6%%points) relative to the leaf area (85%%points), which leads us to the conclusion that the silhouette area estimates are relatively robust against changes in voxel size.

Silhouette to total area estimates
We observed STAR crown ranging between 0.075 and 0.25 (Fig. 6), however, due to the above mentioned limitations in estimating silhouette and total area, we interpreted the results only in relative terms between the individual trees. Overall there was a positive relationship between STAR and voxel size caused by the assumption of random distribution of elements within voxels by the MCF. As explained in Section 4.1, increasing voxel size likely resulted in increased clumping to occur within voxels, and the assumption of random distribution of elements within voxels was thus violated. In addition, small trees (trees 2-5 and 9) exhibited a distinct peak at moderate voxel sizes, similar to what we observed for silhouette areas. It was likely caused by processing artifact introduced by the voxel grid on the small crown volume, as discussed in Section 4.2.
The relative level of STAR crown differed strongly between trees, which is most notable at the smallest voxel size (5 cm). Sparse trees, such as trees 4 and 9, exhibited the highest STAR values and were thus the least clumped trees. On the other hand, trees with crown parts of very high leaf area density (e.g. trees 3, 5, 12, 13) form the opposite extreme, as they were the most clumped. Visually, these two extremes can be observed intuitively from Fig. 1, where the dense and more clumped trees had a distinct volume in the lower part of the crown with an extremely high density relative to the other parts of the crown (cf. the lower third of the crown of tree 3).
Despite the possibility to visually interpret clumping of vegetation elements, it is vital to study their spatial distribution explicitly. That is, without a spatially explicit measure of the crown or stand silhouette to needle silhouette area (i.e. the clumping index 4STAR), a potentially considerable amount of uncertainty and bias is introduced by predicting light interception in canopies using the Lambert-Beer law. LiDAR technologies are currently the only viable option for studying the threedimensional patterns of vegetation elements explicitly. However, LiDAR methods have the above mentioned limitation of the theoretical methods diverging somewhat from the actual measurement process, namely by the effective beam footprint effect, and the lack of a robust method to determine a suitable voxel size for estimating attenuation coefficient from point clouds. Overcoming these limitations would allow the exact analysis of absolute level and variation of STAR for individual tree crowns and entire forest stands. Ultimately, such an improved knowledge about STAR would enable modelers to partly or fully abandon the concept of a Poisson canopy and replace it with more realistic spatial probability distributions, or at the least, justify the use of Poisson canopy for some tree species or stand compositions.

Conclusions
We developed a method to estimate the silhouette to total area ratio from multiple view TLS point clouds, and compared our method with reference measurements in young spruce trees. Our method to estimate STAR at tree scale is based on estimating the spatially explicit attenuation coefficient on a voxel grid, which is further processed to calculate total leaf area, and tree silhouette area in all directions of a sphere. The method we used for estimating the attenuation coefficient showed a tendency to overestimate both leaf and silhouette area relative to reference measurements, with a strong dependence on the spatial scale (i.e. the voxel size). We analyzed our results with regard to the influence of voxel size, and pointed out important error sources that influence the estimation accuracy for attenuation coefficient.
The STAR we observed in young spruce trees shows a significant level of clumping at small scales between 10 and 25 cm, which indicates that the spatial distribution of shoots in a tree crown is not random. In future work, we plan to apply our proposed method to larger trees and forest stands to empirically analyse the degree of clumping in tree crowns and stands. Furthermore, additional research is required to improve the estimation of attenuation coefficient in voxel grids, particularly for the influence of the effective beam size and the consequent lack of resolution of conifer needles.
Finally, a potential bias source we analyzed in this study is clumping. We showed that for large voxel sizes, within-voxel clumping plays an increasingly important role as error source in estimating attenuation coefficient. Therefore, we recommend to keep the voxel size for estimating attenuation coefficient and leaf area as small as possible to avoid the introduction of additional errors that may partly or wholly cancel out other errors inherent to the modified contact frequency or the TLS measurements.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.