Opaque voxel-based tree models for virtual laser scanning in forestry applications

Virtual laser scanning (VLS), the simulation of laser scanning in a computer environment, is a useful tool for field campaign planning, acquisition optimisation, and development and sensitivity analyses of algorithms in various disciplines including forestry research. One key to meaningful VLS is a suitable 3D representation of the objects of interest. For VLS of forests, the way trees are constructed influences both the performance and the realism of the simulations. In this contribution, we analyse how well VLS can reproduce scans of individual trees in a forest. Specifically, we examine how different voxel sizes used to create a virtual forest affect point cloud metrics (e.g., height percentiles) and tree metrics (e.g., tree height and crown base height) derived from simulated point clouds. The level of detail in the voxelisation is dependent on the voxel size, which influences the number of voxel cells of the model. A smaller voxel size (i.e., more voxels) increases the computational cost of laser scanning simulations but allows for more detail in the object representation. We present a method that decouples voxel grid resolution from final voxel cube size by scaling voxels to smaller cubes, whose surface area is proportional to estimated normalised local plant area density. Voxel models are created from terrestrial laser scanning point clouds and then virtually scanned in one airborne and one UAV-borne simulation scenario. Using a comprehensive dataset of spatially overlapping terrestrial, UAV-borne and airborne laser scanning field data, we compare metrics derived from simulated point clouds and from real reference point clouds. Compared to voxel cubes of fixed size with the same base grid size, using scaled voxels greatly improves the agreement of simulated and real point cloud metrics and tree metrics. This can be largely attributed to reduced artificial occlusion effects. The scaled voxels better represent gaps in the canopy, allowing for higher and more realistic crown penetration. Similarly high accuracy in the derived metrics can be achieved using regular fixed-sized voxel models with notably finer resolution, e.g., 0.02 m. But this can pose a computational limitation for running simulations over large forest plots due to the ca. 50 times higher number of filled voxels. We conclude that opaque scaled voxel models enable realistic laser scanning simulations in forests and avoid the high computational cost of small fixed-sized voxels.


Introduction
Metrics derived from airborne laser scanning (ALS) point clouds relate to important structural forest characteristics and can be used to complement and enhance time-consuming conventional forest inventories (FIs) (Maltamo et al., 2014). For predicting forest and tree characteristics from point cloud data, it is important to understand the relation between ALS metrics and FI variables, and the robustness of ALS metrics to different acquisition parameters (e.g., flying altitude, scan angle, pulse repetition frequency, beam divergence). These questions have been addressed by comparing laser scanning data obtained with varying acquisition settings with in situ reference information (Chasmer et al., 2006;Hopkinson, 2007;Morsdorf et al., 2008;Naesset, 2009). Since ALS acquisitions are expensive, only a limited number of configurations can be tested, and usually, the effect of only a single or few acquisition variables is investigated. This makes it difficult to disentangle the multiple and also co-dependent influences of sensor and platform configurations on the point cloud characteristics. To overcome this limitation, extensive sensitivity analyses can be conducted by the simulation of laser scan acquisitions in a virtual environment (Disney et al., 2010;Holmgren et al., 2003;Roberts et al., 2020), hereafter referred to as virtual laser scanning (VLS). VLS allows to simulate practically unlimited amounts of data with different configurations at low cost. The simulated data have perfect ground truth, with pre-defined errors in georeferencing and strip alignment. VLS data may therefore be used to investigate the influence of different sensor and platform configurations on the point cloud metrics. Furthermore, VLS data with perfect ground truth may serve as training and testing data for algorithm development and machine learning in application fields such as tree detection, tree segmentation and forest gap detection.
Essential inputs for laser scanning simulations comprise the sensor and platform configurations and the virtual scenes. To assemble 3D virtual forest scenes, 3D models of trees are needed. Qi et al. (2017) presented a system to reconstruct a large 3D forest scene (1 km × 1 km) by combining a database of representative computer-generated tree models with ALS data of a forest plot. A related approach to generate synthetic forestry data was suggested by Fassnacht et al. (2018) and Schäfer et al. (2019), where tree positions and characteristics are determined using a forest growth simulator and forests are then built by voxelising real 3D point clouds of trees with matching characteristics, stored in a tree database.
Single tree modelling techniques can be categorised by their data source, methodology, and by the geometric primitives used to represent tree components, i.e., the stem and the crown, or individual branches and leaves. In a simple form, crown archetypes such as ellipsoids, hemiellipsoids or cones have been used to represent crown volumes in laser scanning simulations or for modelling canopy reflectance (Calders et al., 2013;Holmgren et al., 2003;Lovell et al., 2005;Ranson and Sun, 2000;Strahler and Jupp, 1990;Widlowski et al., 2014). The crown archetypes can be modelled as solid, i.e., fully opaque (Holmgren et al., 2003), or penetrable, e.g., as turbid medium. Typically, penetrable archetypes are assumed to be filled with foliage which is characterised by structural parameters (leaf area density, angle distribution and size) and spectral properties (reflectance, transmittance) (Calders et al., 2013;North et al., 2010). The choice of these parameters is often based on heuristics.
A further option is to use very detailed synthetic 3D mesh models generated by computer algorithms included in software such as Arbaro (Weber and Penn, 1995), OnyxTREE (Onyx Computing, 2020) or AMAPstudio (Griffon and de Coligny, 2014). Such models have often been applied for virtual laser scanning (Hämmerle et al., 2017;Qi et al., 2017;Qin et al., 2017;Roberts et al., 2020).
Alternatively, trees can be modelled from real 3D point cloud data. Multi-scan terrestrial laser scanning (TLS) point clouds are characterised by high point density and high accuracy, making them a suitable data source for data-driven high-resolution tree modelling. TLS point clouds have been used for direct reconstruction of tree branching structure and distribution, orientation, size, and shape of foliage in sophisticated algorithms (Côté et al., 2011;Hackenberg et al., 2014;Raumonen et al., 2013;Xie et al., 2018). Typically, TLS point clouds acquired under leaf-off conditions are used for the generation of such architectural or quantitative structure models (QSMs). In case of leaf-on data, foliage must be removed in a pre-processing step. This results in lower quality point clouds due to the limitations of leaf removal algorithms, and occlusion in the data caused by leaves, especially towards the top of the trees (Calders et al., 2018). To then obtain leaf-on tree models, leaves or needles have to be added again in a separate step, which requires knowledge of foliage parameters like leaf shape (Milenković et al., 2012), and the distribution of leaf area density, leaf size, and leaf orientation (Åkerblom et al., 2018).
Tree canopies can be modelled from TLS (or other) point clouds in a more straightforward way by dividing the space into a 3D grid of volume elements (voxels). Filled voxels, i.e., voxels with laser returns inside, may then be defined as solid or turbid medium and assigned spectral properties (Gastellu-Etchegorry et al., 2015;Widlowski et al., 2014). Having defined a threshold for the minimum number of points for a voxel to be considered as filled, tree reconstruction using regular opaque cubic voxels requires only one more parameter, the voxel size. This approach can be used with both leaf-off and leaf-on data. Due to its simplicity, voxelisation is useful when exact branching structure and foliage representation are not required and large forest areas have to be modelled.
The choice of the tree modelling approach influences the simulated point cloud, the accuracy of the derived metrics, and the computational cost of the simulation. On the one hand, very detailed tree reconstruction with many small primitives, e.g., triangles or voxels, can be limiting for the simulation of large forest areas due to memory and runtime constraints. On the other hand, oversimplified tree models lead to less realistic point clouds. This motivates our in-depth investigation of the influences of different opaque voxel tree representations on the structural properties of simulated laser scanning point clouds.
Our objectives are (i) examining the suitability of opaque voxelbased reconstruction methods from TLS point clouds for realistic airborne and UAV-borne laser scanning simulations over forests and (ii) designing a voxel-based tree model that aims to decouple the number of opaque voxels from the level of detail of the model using a scaling procedure. The overarching aim is to develop an approach that allows simulation of laser scanning point clouds which are geometrically consistent with real point clouds at low computational cost.
We compare three approaches to create opaque voxel-based forest representations: (1) regular voxels examined with different voxel sizes (2) dynamically scaled voxel cubes using voxel plant area density as scaling reference, which is estimated from the TLS point cloud (3) dynamically scaled voxel cubes with additional shifting of the cubes towards the centre of gravity (CoG) of TLS points within the initial voxel space.
VLS is performed over the different voxel-based forest representations in one airborne and one UAV-borne scenario. Modelling and simulation accuracy are evaluated by comparing simulated point cloud metrics and tree metrics of individual trees to the metrics of corresponding point clouds obtained with similar configurations in real-world acquisitions. Our first hypothesis is that the scaling (Model 2) leads to a detailed representation of the crown without the need of a detailed voxel grid, i. e., a high number of small voxels. Secondly, we hypothesise that shifting the scaled cubes towards the CoG (Model 3) further improves the accuracy of the point cloud metrics and tree metrics because the resulting representation more closely resembles the actual tree structure in terms of location of potential scatterers, and any influence from the regular grid pattern is reduced.
In the following, we present our study area and datasets of real-world acquisitions (Section 2). Subsequently, we explain our methods for forest modelling, laser scanning simulation and metric accuracy analysis (Section 3). We present our results in Section 4 and discuss our findings in Section 5, before drawing conclusions (Section 6).

Study area and datasets
This study utilises TLS point clouds for forest reconstruction and subsequent VLS. VLS point cloud metrics and tree metrics are then compared to metrics derived from airborne laser scanning (ALS) and UAV-borne laser scanning (ULS) point clouds which were obtained with similar configurations in the field as reference dataset.

Study area
The study area is the communal forest of Bretten in Baden-Württemberg, Germany (49 • 0 ′ N, 8 • 42 ′ E), in which groups of trees were surveyed with TLS in six small rectangular plots (Fig. 1). The plots cover areas between 0.19 ha and 0.5 ha. 30 selected single trees within the plots, which represent the range of tree characteristics of the forest stand well, are investigated in the detailed validation procedure (Section 3.4). These trees are referred to as our target trees. Half of the trees are broadleaf deciduous trees (13 Fagus sylvatica, 2 Carpinus betulus), the other half are coniferous trees (15 Picea abies). The ranges, means and medians of tree height, crown base height (CBH) and crown projection area (CPA) are presented in Table 1. CBH and CPA were derived from TLS point clouds and tree height from ULS point clouds as described in Section 3.4.

TLS data acquisition and pre-processing
TLS point clouds were acquired with a RIEGL VZ-400 scanner mounted on a tripod under leaf-on conditions in June and July 2019. The laser scanner operates at a wavelength of 1550 nm and with a beam divergence of 0.3 mrad (measured at the 1/e 2 points, RIEGL Laser Measurement Systems, 2017). Scans were performed with a pulse repetition frequency (PRF) of 300 kHz and vertical and horizontal angle increments of 0.017 • (0.3 mrad), resulting in 3 mm point spacing at 10 m range for a single scan. In each of the six campaigns, a group of trees was surveyed from five to eight scan positions, distributed on a circle or ellipse around those trees (Fig. 1). Each scan had a horizontal field of view between 60 • and 100 • and a vertical field of view of 100 • . Where the trees of interest were very high in relation to their distance to the scanner, vertically tilted scans were added to the scan position to increase the vertical coverage of the point cloud. Five cylindrical reflective targets mounted on tripods were placed around the trees in the centre of the scene. Further circular reflectors were pinned on the target trees to help with the tree identification in the point cloud and with the co-registration of the single scans. For each TLS campaign, RTK GNSS measurements from one scan position, referred to as the main stable position, and one tiepoint were used to georeference the respective scan. Remaining scans were co-registered to the main position using reflective targets. This tiepoint-based registration was carried out in RiSCAN PRO (Version 2.8.0). A fine alignment of scans was performed with the Multi Station Adjustment function in RiSCAN PRO, which uses a variant of the iterative closest point (ICP) algorithm based on planar areas in the point cloud. Lastly, a co-registration to the ULS dataset was performed to  accurately georeference the TLS data. This was necessary because of the low accuracy of the in situ GNSS measurements due to the weak signal in the forest. For this registration step, tree stem positions were computed from the downsampled TLS and the ULS point clouds by calculating local normal vectors, linearity and 3D point densities using OPALS (Version 2.3.1, Pfeifer et al., 2014). The TLS-derived stem positions were then manually matched to the ULS-derived stem positions in 2D and the resulting rigid transformation was applied to the TLS point cloud. After estimating and correcting the height offset by selecting corresponding points in a vertical profile in both point clouds, a fine alignment was applied with the ICP implementation in OPALS.

ALS and ULS data acquisition
Airborne laser scanning was conducted over the entire Bretten forest area on 5 July 2019 with a RIEGL VQ-780i mounted on an aircraft of type Cessna C207. The scanner has a beam divergence of 0.25 mrad (measured at the 1/e 2 points, RIEGL Laser Measurement Systems, 2019) and was scanning with a PRF of 1000 kHz and a scan frequency of 225 lines per second at approximately 650 m altitude above ground level. The scan angle was ±30 • and the flight strip spacing of 175 m resulted in about 76% strip overlap on the ground. The airplane was flying with approximately 100 knots (51 m s − 1 ). The resulting point cloud has a pulse density (i.e., density of last returns) of 70 to 80 points/m 2 .
UAV-borne laser scanning point clouds were collected with a RIEGL miniVUX-1UAV mounted on a DJI Matrice 600 Pro on 24 August and on 12 September 2019 (DJI, 2018; RIEGL Laser Measurement Systems, 2020). The UAV was moving with around 5 m s − 1 at 60 m to 80 m altitude above ground level. The scanner has a beam divergence of 1.6 mrad × 0.5 mrad (measured at 50% peak intensity, RIEGL Laser Measurement Systems, 2020). This means the scanner has an elliptical footprint in the form of a two-dimensional Gaussian, where the energy drops to 50% of the maximum at 0.8 mrad off the centre for the semi-major and at 0.25 mrad off the centre for the semi-minor axis. PRF was fixed at 100 kHz and scan frequency was set to 50 lines per second. The scan angle was ±90 • . The acquisition was performed with two overlapping double grids, where the spacing between the scan lines was 26 m to 30 m (Fig. 5). The resulting point cloud has a pulse density of 800 to 1200 points/m 2 .

Methods
The main steps in this study comprise tree segmentation, tree modelling, laser scanning simulation and tree metric computation, followed up by an accuracy assessment. The workflow is visualised in Fig. 2.

Tree segmentation
Individual trees are automatically segmented with the CompuTree software (Version 5.0.054b, Computree Group, 2017) using a Euclidean clustering approach combined with a competitive Dijkstra's algorithm. The methods are part of the SimpleTree Plugin (Version 4.33.06 Beta, Hackenberg et al., 2015;Hackenberg, 2017). To obtain accurate 3D tree delineations, the automatic segmentation of target trees is improved manually. This is done by visual assessment and manual digitisation. The automatically extracted point clouds of target trees are compared with neighbouring point clouds. Wrong points are removed from the target tree point cloud and points wrongly assigned to neighbouring tree point clouds are added. We perform this manual correction iteratively by two editors to ensure high quality tree delineation. Extracted tree point clouds are filtered to keep only points with a pulse shape deviation (as defined by RIEGL) lower than 50, i.e., returns with extremely high echo width are removed. For our accuracy analysis of single tree metrics, the point clouds of 15 coniferous trees and 15 deciduous broadleaf trees, our target trees, are selected.
The segmented TLS point clouds of the target trees serve as templates to extract point clouds of the same trees from the reference ALS and ULS datasets. For each TLS point in the tree point cloud, all points in a neighbourhood of 0.5 m (ALS) and 0.3 m (ULS) are queried from the ALS and ULS datasets and labelled as belonging to the tree. Errors resulting from this automatic procedure, i.e., labelling of points that obviously do not belong to the respective tree, are manually corrected.
The full TLS point clouds are clipped to smaller rectangular plots containing the target trees ( Fig. 1). To adequately capture occlusion effects by surrounding trees in the laser scanning simulations, circular buffers of 20 m are created around each target tree position. The extent of the six rectangular forest plots is determined as the minimum axisaligned bounding rectangles enclosing the 20 m tree buffers. The

Voxel-based tree and forest modelling
To create virtual scenes for laser scanning simulations, the TLS point clouds of the target tree stems are meshed and the TLS point clouds of the target tree crowns and all surrounding trees are voxelised with the methods summarised in Table 2 (Fig. 3). Our reasoning for using a hybrid approach for the target trees (i.e., crowns are voxelised while stems are meshed) is that with a high resolution multi-scan TLS setup, small point spacing and low occlusion effects allow for effective reconstruction of the tree stems. Furthermore, tree stems are assumed solid and unable to be penetrated by the laser beam, unlike the tree canopy where gaps are present. A Poisson surface reconstruction (Kazhdan et al., 2006;Kazhdan and Hoppe, 2013) implemented in CloudCompare (Version 2.10.2) is performed for the stem point cloud, resulting in a stem mesh.
We compute simple fixed-sized crown and forest models with four different voxel sizes (i.e., voxel side lengths): 0.25 m, 0.1 m, 0.05 m and 0.02 m ( Table 2). Each cubic voxel cell containing at least 1 point is considered a filled voxel providing a surface to be scanned in the laser scanning simulations. Both leaves and branches may be present within a voxel.
For the scaled voxel models, point clouds are voxelised with the AMAPVox software (Version 1.5.5, Verley et al., 2019). The co-registered RiSCAN PRO projects are given as input to AMAPVox to retrieve the pathway of each emitted laser pulse. Transmittance, or gap probability P gap , is then computed using Eq. (1) (Vincent et al., 2021): where P gap is the mean weighted gap probability, BFOut i the outgoing beam surface fraction of pulse i, BFEnt i the entering beam surface fraction of pulse i, S i the cross section of the pulse at the voxel centre, and l i the potential length of the optical path of pulse i. The beam fractions are weighted according to the echo rank: The residual fraction of energy loss of a pulse is set to the inverse of the number of targets (i.e., 50% in case of two targets, 33% in case of three targets, etc., Vincent et al., 2017). Following Beer-Lambert law, local plant area density (PAD) can be calculated from transmittance (Eq. (2), Grau et al., 2017, Béland et al., 2014b, Vincent et al., 2017.
where G(θ) is the leaf projection function (Nilson, 1971), which can be computed for given leaf angle distributions. The term plant area density is used here instead of leaf area density, defined as the one-sided leaf area per unit volume (m 2 /m 3 ), because we do not discriminate between leaves and wood (Vincent et al., 2017). For the computation of G(θ), the theoretical one-parameter distribution functions by de Wit (1965) are used to describe the distributions of leaf angles within the crowns (Wang et al., 2007). Following the results of Liu et al. (2019), Wagner and Hagemeier (2006) and Zhu et al. (2018), a planophile distribution is assumed for the broadleaf trees (F. sylvatica and C. betulus). For the conifers (P. abies), an erectophile distribution is assumed according to Zhu et al. (2018). The voxel resolution is set to 0.25 m to make sure that multiple leaves and twigs are present in a single voxel, i.e., the voxel size is larger than the leaves (Schneider et al., 2019;Béland et al., 2014a). All estimates larger than a defined PAD max , are limited to PAD max to exclude unlikely high values as outliers. From these estimates, we scale each voxel into a cube by its PAD value using Eq. (3): where a is the cube side length in m, a 0 the base voxel size (i.e., grid size) in m and α the scaling factor, here set to 0.5, ensuring that the surface area (A = 6a 2 ) of the scaled cubes is proportional to estimated PAD. PAD max is set to 5.0 m 2 /m 3 for deciduous trees and 8.0 m 2 /m 3 for coniferous trees. This way, voxels of coniferous trees are scaled to smaller cubes, accounting for the small needles in the crown, making it more transmittive. The values were chosen empirically by comparing simulation outputs for different values for PAD max . This variably sized voxel model is used to examine our hypothesis that scaling leads to a detailed representation of the crown without the need for a detailed voxel grid, i.e., a high number of small voxels (cf. Section 1).
For the investigation of Hypothesis 2, we shift the scaled cubes from their original positions at the voxel cell centres to avoid an artificial regular structure and to resemble the input TLS data (i.e., the location of the scatterers) more closely. For this, we calculate the CoG of all TLS points located in the respective voxel cell and move the centre of the scaled cube towards the CoG. Shifting is only performed within the voxel cell, so scaled voxels are not allowed to overlap. To investigate the effect of the shift on the resulting point cloud, a random shift of cubes within the voxel space is applied for comparison. This is achieved by multiplying each component of the cubes' maximum possible shift vector (i. e., not shifting it outside of the voxel space) with an independent random number drawn from a standard uniform distribution.
Exemplary tree models of different species reconstructed with 0.25 m fixed-sized voxels, and scaled and shifted voxels are shown in Fig. 4 together with the associated TLS point clouds.

Laser scanning simulations
Laser scanning simulation is conducted with the Heidelberg LiDAR Operations Simulator HELIOS++ (Version 1.0.8, Winiwarter et al., 2021a,b). The six scenes consist of the reconstructed rectangular TLS forest plots with the target trees located in the respective centres (cf. Figs. 1 and 3). In both simulations with scaled and fixed-sized voxels, voxels are opaque, which means simulated beams are reflected from the outside of one of the faces of a cubic voxel. The ground is represented by a simple planar surface at the approximate elevation of the terrain. In the laser scanning simulations, each simulated echo is attributed to the object from which it was returned by an ID. By modelling each target tree stem and crown as individual objects, we can easily extract the simulated target tree point clouds by these object IDs.
Scenes are surveyed with HELIOS++ in one ULS and one ALS scenario with the settings shown in Table 3. The configurations are based on the real laser scanning campaigns which are used as reference in the accuracy assessment. In the ALS scenario, flying altitude above ground level is 658 m to 725 m. In the ULS scenario, flying altitude above ground level is 64 m to 72 m. Exemplary flight trajectories for P3 are visualised in Fig. 5. In the ULS simulations over P2, one of the two overlaid double cross patterns is flown at 4 m s − 1 instead of 5 m s − 1 . HELIOS++ does not support the simulation of the elliptical beam footprint of the RIEGL miniVUX-1UAV used in the reference campaigns. In the simulations, we therefore assume that the scanner has a circular beam footprint. Beam divergence is set to the divergence of the minor axis of the elliptical beam of the real scanner. The distance between the flight lines is around 175 m in the ALS scenario and around 26 m to 30 m in the ULS scenario. We investigate the number of first returns in the target tree point clouds to make sure the simulated platforms and beam deflectors work as expected and produce the same number of pulses as in the real reference surveys (Appendix A, Figure A1).
The simulation framework HELIOS++ approximates beam divergence by performing ray tracing on multiple subrays within the laser cone of given divergence (Winiwarter et al., 2021b). The number of subrays is controlled with the "beam sample quality" parameter, defining the number of concentric circles on which subrays are sampled. For the simulations in our study, beam divergence is simulated in three circles resulting in 19 subrays. In the simulated waveform, returns are detected as local maxima in a moving temporal window given by the "window size" parameter in nanoseconds (Winiwarter et al., 2021a). A local maximum will be detected only if there is no larger value within [− window size; +window size] from the current position. By multiplying by the speed of light, we can therefore convert the size of the time window to the vertical dead zone (minimum distance between two detected points within a beam) in metres. We use a window size of 1.75 ns for the ALS scenario and of 1.5 ns for the ULS scenario, corresponding to vertical dead zones of 0.525 m and 0.45 m. The window sizes are chosen to match the multiple return ratios observed in the real laser scanning tree point clouds. Smaller window sizes lead to a higher number of detected returns, while larger window sizes result in a lower   number of detected returns. The z-coordinates of extracted simulated point clouds of the target trees and reference point clouds from real laser scanning are normalised with the same value for matching trees. This way, heights of simulated and real tree point clouds can be directly compared.

Accuracy analysis of point cloud metrics and tree metrics
Selected point cloud metrics and tree metrics are computed from the VLS point clouds and the real reference point clouds of the target trees (Table 4 and Fig. 6). These metrics are compared between the simulated and the reference point clouds to derive the accuracy that can be achieved with the different voxel modelling approaches.
The relative frequency of multiple returns gives information about laser penetration into the tree canopies (Takahashi et al., 2006), especially in combination with the mean height H mean , the standard deviation of height H sd , and the minimum height H min . Here, we define H min as the height of the lowest return on the respective target tree. Height metrics (H mean , H sd , H min and height percentiles) and density metrics are affected by tree characteristics (tree height, crown height, crown area) and foliage density (in real acquisitions) or density of the voxel canopy representation (in simulations).
In addition to these statistical point cloud metrics, we derive tree metrics for the single target trees (Table 4 and Fig. 6). Tree height is obtained as the maximum height above ground for each tree point cloud.
To determine the crown base height (CBH), the point cloud is divided into height sections of 0.1 m. In each section, the maximum horizontal distance between any two points is calculated. Starting from the bottom, we iterate over the height sections until the maximum horizontal point distance exceeds a threshold of 1.5 m. The height of the centre of this section is defined as the CBH (Fig. 6). Due to the fixed section height, the algorithm as implemented only achieves decimetre resolution. It detects the lowest branch exceeding a certain length, but does not consider if it is connected to the rest of the crown. CBH is an important parameter for forest fire simulations (Gajardo et al., 2014). It is also the basis for crown-related metrics like crown projection area, crown depth, and crown volume. Crown projection area (CPA) is calculated as the area of the planar (x-y) convex hull of all crown points, i.e., all points above the CBH. All metrics are computed on the single tree point clouds only, excluding ground points and points recorded on the surrounding trees. We use all returns, i.e., first and all subsequent returns per pulse.

Results
In this section, we present the agreement of metrics derived from simulated point clouds with those derived from real (reference) point clouds of the same trees, i.e., our target trees, for the different voxel modelling approaches. Where needed, a differentiation is made not only between the flight scenarios, but also between the tree species. A complete collection of metrics for each tree derived from the different scenarios can be found in Appendix A, Detailed listing of metrics for each single tree.

Height and density metrics
When using models with the larger fixed-sized voxels, i.e., 0.1 m to 0.25 m, mean height percentiles are higher than reference mean height percentiles (Fig. 7). For these two largest voxel sizes, mean point densities are overestimated in the upper two height bins and underestimated in height bins 1 to 8 (Fig. 8). Correspondingly, larger voxels lead to overestimation of H mean (Appendix A, Figure A3) and H min , and underestimation of H sd (Fig. 9), because most simulated pulses are reflected at the top of the voxel canopy. Models with smaller voxels result in lower height percentile curves, more returns in lower height bins and lower H min . For broadleaf trees, the mean point density in the upper height bin is underestimated with the smallest voxel size, whereas point densities in height bins 8 and 9 are overestimated. For coniferous trees in the ALS scenario, all four voxel sizes result in overestimation of the average height percentile curves and the mean point density in the upper two height bins ( Fig. 7b and Fig. 8b, blue lines/bars).
For all voxel models, overestimation of H min is more pronounced in the ALS scenario compared to the ULS scenario due to the lower pulse density of ALS acquisitions. For broadleaf trees, there is one outlier, where H min is between 0.5 m and 2.5 m when simulated with the two smaller fixed-sized and the scaled voxels but at around 18 m in the reference data ( Fig. 9c and Appendix A, Figure A4). This demonstrates that also in the reference point clouds, the lowest tree return may not be very close to the ground. Underestimation of this magnitude does not occur in the ULS scenarios because the trunk base of the tree is captured well in all reference point clouds.
With scaled voxels, height percentile curves are close to the reference for broadleaf trees but, like with all other voxel models, overestimated for coniferous trees (Fig. 7). H sd is underestimated with scaled voxels by around 10%, but absolute errors of H sd are lower for scaled voxels than for all fixed-sized voxel models (Fig. 9). Errors of H min for scaled voxels are very small and comparable to the smallest 0.02 m fixed-sized voxels. Shifting scaled cubes randomly has a negligible effect on height and density metrics. Shifting to the CoG slightly decreases the height percentiles and the mean height (Appendix A, Figure A2).

Multiple returns
Compared to crown models made up of large fixed-sized voxels, models made up of small fixed-sized or scaled voxels allow for a higher frequency of multiple returns, i.e., higher crown penetration, because these models have more gaps. With the smallest 0.02 m voxels, the frequency of multiple returns is overestimated in the ULS scenario by around 8 and 20 percentage points for coniferous trees and broadleaf trees, respectively (Fig. 10). In the ALS scenario, multiple returns are overestimated for broadleaf trees but underestimated for coniferous trees with the 0.02 m voxels. With scaled voxels, the frequency of multiple returns is overestimated for broadleaf trees in the ALS scenario and coniferous trees in the ULS scenario.
In both ALS and ULS reference point clouds, the frequency of multiple returns is significantly higher for coniferous trees (ALS: M = 73.4%, ULS: M = 33.1%) than for broadleaf trees (ALS: M = 55.9%, ULS: M = 28.1%), t(28) = 5.54, p = .000006 for ALS and t(28) = 2.31, p = .029 for ULS. Simulations with fixed-sized voxels do not reflect these interspecific difference. Simulations with scaled voxels however consistently yield a higher number of multiple returns for coniferous trees than for broadleaf trees, which is in agreement with the reference. While random shifts do not influence the frequency of multiple returns, shifting the voxels towards the CoGs reduces the frequency of multiple returns in both scenarios and for both species.

Tree metrics
Overestimation of CBH, and underestimation of CPA increase with larger voxel size (Fig. 11). Because crowns reconstructed with large fixed-sized voxels are less penetrable, they introduce artificial occlusion effects. Most simulated laser pulses are reflected at the top of the voxel canopy and few pulses reach the lower crown or the stems. CBH estimates agree best with the reference values for the more penetrable 0.02 m fixed-sized voxels and the scaled voxels. CPA is estimated most accurately with scaled voxels (Fig. 11, Appendix A, Figures A6 and A7). CPA errors increase if the scaled voxel cubes are shifted to the CoGs.
Tree height errors are small across all employed voxel models (Fig. 12). For fixed-sized voxels, a reduction of the median estimated tree height can be observed with decreasing voxel size (Fig. 12). In the ALS scenario, there is an outlier amongst coniferous trees at different voxel sizes, for which tree height is underestimated by more than 0.8 m.
For coniferous trees, tree height is overestimated more with scaled voxels than with the smallest fixed-sized voxels. Some of the tree height errors can be traced back to the differences between TLS-, and ALS-or ULS-derived tree heights (Fig. 12b and c). Because TLS point clouds were used to create the models, these differences influence the tree height of simulated point clouds.

The influence of voxel size on occlusion and derived tree metrics
The size of the opaque voxels used for tree modelling has a major influence on the frequency of gaps, which is consequently influencing the vertical point distribution in the simulated point clouds. For larger voxel sizes, artificial occlusion effects are introduced, which lead to decreasing quality of CBH estimation. Not only the stem, but also parts of the crown may be occluded when using larger voxel sizes. If simulations were performed with isolated treesas opposed to trees embedded in a forest stand like in this study -CPA would tend to be overestimated more significantly with increasing voxel size. However, we observe stronger CPA underestimation with increasing voxel size, which is a result of occlusion effects by overtopping neighbouring trees. Occlusion effects are also present in the real ALS and ULS data, but are weaker because the real tree canopies are more transmittive than the canopies reconstructed with large voxels. The influence of voxel size on CBH and CPA estimates is illustrated in Fig. 13 by comparing ALS point clouds generated in simulations employing scaled and 0.25 m fixedsized voxels. Unless using sufficiently small voxel sizes (e.g. 0.02 m for our data), fixed-sized voxels for forest reconstruction may results in little to no returns from tree stems, suppressed trees (i.e., trees with no exposed crown area; Magnussen et al., 1999) and large parts of the ground or understorey.
The lack of points sampled on tree stems in ALS scenarios is not only a result of limited crown penetration due to artificial occlusion, but also of lower pulse density and narrow scan angles compared to ULS, where stems are sampled better in both real and simulated acquisitions.
For one coniferous tree, tree height is underestimated by several decimetres in the ALS scenario (Fig. 12). For larger voxel sizes, this may be explained by occlusion of the tree top by surrounding trees. For very small voxel sizes and scaled (and shifted) voxels, it is more likely that tree height is underestimated because the tree top is missed by the simulated laser beam and/or the area within the beam is too small to create a return.
The relation of increasing accuracy in the calculated metrics with decreasing voxel sizes is limited to a certain voxel threshold. If voxel sizes are below this threshold, the model will become too transparent and will have too many gaps which the laser beam can pass through. This especially influences the point cloud metrics (height metrics and percentiles, densities, multiple returns). A lower limit for suitable voxel sizes may be determined from the point spacing in the TLS point cloud.
To some extent, too small voxels are also a problem in our study, where with 0.02 m fixed-sized voxels the frequency of multiple returns is overestimated substantially in the ULS simulation. Furthermore, crown penetration into the top of the canopy is overestimated for broadleaf trees, which is reflected by the much lower point densities in the upper height bin for 0.02 m fixed-sized voxels (Fig. 8). This can be explained by characteristics of the TLS input point clouds, which we discuss in the next section.

Limitations of TLS input point clouds
Because of the scanning geometry from below the canopy, TLS point clouds have lower point density towards the top of the canopy due to occlusion effects caused by leaves (Wilkes et al., 2017). This is also reflected in TLS-derived tree height, which is underestimated compared to tree height from ULS reference point clouds for broadleaf trees (Fig. 12c).
The problem of incomplete sampling of the top of the canopy will be reflected in the tree models (Fig. 4c) and in the simulated point clouds. Reconstructed with very small or scaled voxels, the model will be more transparent in the top of the canopy than in the lower canopy. When simulating a ULS flight over a broadleaf tree model created from TLS, the virtual ULS data will show lower tree height and higher penetration into the top of the canopy than real ULS data at the same location when using a very detailed modelling approach (Fig. 12c). This effect is species dependent and should be kept in mind for VLS-based studies. For coniferous trees, TLS tree height estimates are similar to reference ULS Fig. 10. Mean frequency of multiple returns of target tree point clouds of broadleaf trees ("B") and coniferous trees ("C") simulated with different tree representations in the two different scenarios. Black bars and lines = Reference. Blue bars = Simulated with fixed-sized voxels with different side lengths. Red bars = Simulated with scaled voxels with no shifts, shifts to CoG or random shifts. Note that the simulation was optimised to achieve reliable estimates of the frequencies of multiple returns by adapting the "window size" parameter (Section 3.3). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) estimates and higher than ALS estimates. TLS point clouds of coniferous species are less affected by occlusion than those of deciduous trees under leaf-on conditions. This is due to the differences in (a) foliage structure with conifer canopies being more penetrable than broadleaf canopies and (b) the crown shape as coniferous trees have a more narrow crown than deciduous trees in our study area.

Tree height estimates
Comparing the different real datasets of our target trees (TLS, ULS, ALS), average tree height estimated from ULS is the highest. This results from the combination of high sampling density and small beam footprint compared to airborne acquisition, and the scanning geometry from above. We hence consider ULS-derived tree heights as most reliable. Average tree height from ALS is underestimated for all tree species. In airborne acquisitions, the area of the tree tops illuminated by the laser beam may be too small to create a return signal strong enough to be detected and the beam penetrates further into the canopy (Andersen et al., 2006;Gaveau and Hill, 2003;Lefsky et al., 2002;Wieser et al., 2016).
We reproduced this tendency for ALS tree height underestimation in our simulations for coniferous trees: Tree heights from ALS data simulated over fixed-sized voxel models with fine voxel grids (0.05 m and 0.02 m) are very close to the ALS reference and lower than TLS tree heights, which were the basis for the tree modelling (Fig. 12b). This is not the case for downscaled voxels, presumably because of relatively large (i.e., dense) voxels near the tree tops (voxel sizes >0.05 m).

Complexity of model generation vs. performance of simulations
Our obtained metric accuracies suggest that small voxel sizes, i.e., 0.05 m or 0.02 m for our data, are best for realistic simulations when using voxels of fixed uniform size. But the high number of voxels required to reconstruct a forest scene with a fine voxel grid can be a strong limitation for laser scanning simulations in terms of hardware requirements and the computational cost for building the virtual scene.
We suggest downscaling of voxel cubes generated at lower resolution for avoiding this problem. In our approach, we build voxels at a regular grid of 0.25 m and then scale them to smaller cubes using the estimated local plant area density. The number of filled voxels in the resulting models is at least 50 times lower than in the 0.02 m resolution fixedsized model (up to a factor of 100 or more for single target trees) but they achieve equally or more accurate estimates of simulated tree and height metrics. One drawback is that the models are more complex in their construction because they require the computation of PAD. Our approach is dependent on the quality of PAD estimations and on the scaling formula (PAD max and scaling factor α). The trade-offs between modelling complexity, simulation performance (using the proxy "number of voxels") and metric accuracy are summarised in Table 5.

The influence of shifting scaled voxel cubes
The second hypothesis in our study is that additional shifting of the single voxels to the CoGs of TLS points within the original voxel space improves the accuracy of metrics derived from simulated point clouds. The shift reduces the cloud-to-cloud distances between the simulated point clouds and the TLS reference point clouds (Appendix A, Figures A8  to A10). However, only few of the calculated metrics are affected by the shift. The number of multiple returns decreases when shifting voxels to the CoGs. This may be explained by the sorting and clumping of voxels through the shift, reducing the number of partial hits. Furthermore, estimated tree height and CPA decrease, which improves tree height estimates, but deteriorates CPA estimates. Performing a random shift within the original voxel space has negligible influences on all metrics.

Conceptual considerations of the VLS-based approach
Relevant structural tree metrics from virtual ALS and ULS point clouds are predicted with high accuracy when small voxel sizes (i.e., 0.02 m) or scaled voxels are used. In contrast, simulating laser scanning height distribution and multiple returns in tree canopies realistically proves more difficult. For instance, simulated crown penetration using the smallest 0.02 m fixed-sized voxels and the scaled voxels is overestimated for broadleaf trees in both ALS and ULS scenarios but underestimated for coniferous trees in the ALS scenario, shown by the mean point densities in the 10 equal height bins (Fig. 8). This might be related to the quality of the TLS point cloud used for modelling, the modelling approach itself, or the simulation of the beam interaction with the scene. Differences between broadleaf and coniferous trees concerning foliage characteristics, branching structure and crown shape result in different interactions between the laser beam and the trees, as seen, e.g., in the density metrics (Fig. 8) and the frequency of multiple returns (Fig. 10). With voxels scaled by PAD, these differences are better accounted for than with fixed-sized voxels, especially after fine tuning PAD max . The scaling approach depends on the quality of the local TLSestimated PAD values. The problem with the Beer-Lambert law-based (also called gap fraction) method is that increased uncertainty in transmittance rates will lead to overestimation of PAD (Vincent et al., 2021). High sampling variance due to poorly sampled voxels will lead to large uncertainty in local transmittance values. Because PAD is proportional to the logarithm of transmittance, these uncertainties will propagate to PAD. Thus, uncertainties of PAD will be larger for smaller values of transmittance (Vincent et al., 2017). Better results may be achieved with a bias-corrected Maximum Likelihood Estimator as shown by Pimont et al. (2018).
Merging TLS and ULS point clouds for tree and forest reconstruction could potentially alleviate the problems of occlusion of the top of the canopy, and tree height underestimation which are often present in TLS point clouds. However, fusion of TLS and ULS would also increase the demands and efforts for data capturing and processing (e.g., coregistration) significantly.
One inherent property of simulation studies is that influences of the simulation as implemented cannot be separated from influences of the choice of the (forest) scene model. Simulation results are affected by both. For example, the frequency of multiple returns is highly influenced not only by the employed tree model but also by the time window used to extract the maxima from the simulated full waveform.
In our VLS-based concept, we show that the voxel modelling approach has an important influence on simulation output and performance. Our simulated acquisitions over forest scenes reconstructed from TLS with small fixed-sized or scaled voxel models produce realistic VLS point clouds. With the limitations in mind, e.g., concerning the quality of the input point clouds or the need to account for differences between species, such simulations can be used for survey planning, benchmarking and comparison of algorithms, and sensitivity analyses. In forestry, VLS can be of great value for developing algorithms for applications such as single tree detection, single tree delineation or canopy openness estimation. Our study makes an important contribution for establishing VLS as a complementary research and data acquisition tool by validating the simulation output with real data.

Conclusion
This study investigates different opaque voxel-based tree reconstruction methods from TLS point clouds regarding their ability to create models for realistic simulation of ALS and ULS point clouds. This is done by comparing structural point cloud metrics and tree metrics of simulated and real laser scanning data of selected target trees within virtual reconstructed forest plots. In addition to investigating binary voxel models with fixed voxel sizes, an approach of scaling voxels to smaller, variably sized cubes proportional to their plant area density is presented. This enables to voxelise TLS point clouds using a comparably large voxel grid of 0.25 m (i.e., a low number of voxels), yet allowing for a realistic simulation of crown penetration. This is not possible using fixed-sized voxel models at 0.25 m resolution. Such large fixed-sized voxels lead to occlusion effects in the simulated point cloud, causing most returns to occur in the upper canopy and preventing reliable determination of crown base height and crown projection area. The magnitudes of errors for the investigated point cloud metrics and tree metrics are the lowest with scaled voxels or with the smallest fixed-sized voxels (0.02 m). The smallest fixed-sized voxel models however result in approximately 50 times more geometric primitives. Our newly suggested voxel scaling method therefore reduces memory and runtime requirements compared to fixed-sized voxel models with high resolution, i.e., large number of small voxels. Shifting the scaled voxels within the initial voxel volume either towards the centre of gravity of TLS points within each voxel or  randomly does not improve the quality of the metrics. Our results show that the developed scaled voxel approach may be suitable for large-area laser scanning simulations of forest stands due to a balanced trade-off between realism in results and computational cost of simulations.
Such simulations can open up new pathways to fully understand the relationship between forest structure, point cloud data and laser scanning acquisition settings (including flight and scan position planning). This is especially important to further optimise automated operational forest inventories based on airborne and terrestrial laser scanning.

Data statement
The simulated and real point clouds, the files to reproduce the simulations, and the code for preparing the data and deriving the metrics are available through heiDATA, the data repository of Heidelberg University (https://doi.org/10.11588/data/MZBO7T). The point clouds acquired in the field are also part of our data publication with PANGAEA (https://doi.org/10.1594/PANGAEA.933426), albeit with a different naming scheme and in a global coordinate system. For more information, please contact the corresponding author.

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.