<figure><img src="image_1.png" /></figure>

Article

# Influence of 3D Spruce Tree Representation on Accuracy of Airborne and Satellite Forest Reflectance Simulated in DART

Růžena Janoutová <fnref n="1" />,*,† D, Lucie Homolová ¹D, Zbyněk Malenovský ²D, Jan Hanuš ¹D, Nicolas Lauret <fnref n="3" /> and Jean-Philippe Gastellu-Etchegorry ³D

Received: 22 February 2019; Accepted: 19 March 2019; Published: 26 March 2019

<figure><img src="image_2.png" /></figure>

check for
updates

**Abstract:** Advances in high-performance computer resources and exploitation of high-density terrestrial laser scanning (TLS) data allow for reconstruction of close-to-reality 3D forest scenes for use in canopy radiative transfer models. Consequently, our main objectives were (i) to reconstruct 3D representation of Norway spruce (Picea abies) trees by deriving distribution of woody and foliage elements from TLS and field structure data and (ii) to use the reconstructed 3D spruce representations for evaluation of the effects of canopy structure on forest reflectance simulated in the Discrete Anisotropic Radiative Transfer (DART) model. Data for this study were combined from two spruce research sites located in the mountainous areas of the Czech Republic. The canopy structure effects on simulated top-of-canopy reflectance were evaluated for four scenarios (10 × 10 m scenes with 10 trees), ranging from geometrically simple to highly detailed architectures. First scenario **A** used predefined simple tree crown shapes filled with a turbid medium with simplified trunks and branches. Other three scenarios used the reconstructed 3D spruce representations with **B** detailed needle shoots transformed into turbid medium, **C** with simplified shoots retained as facets, and **D** with detailed needle shoots retained as facets. For the first time, we demonstrated the capability of the DART model to simulate reflectance of complex coniferous forest scenes up to the level of a single needle (scenario **D**). Simulated bidirectional reflectance factors extracted for each scenario were compared with actual airborne hyperspectral and space-borne Sentinel-2 MSI reflectance data. Scenario **A** yielded the largest differences from the remote sensing observations, mainly in the visible and NIR regions, whereas scenarios **B**, **C**, and **D** produced similar results revealing a good agreement with the remote sensing data. When judging the computational requirements for reflectance simulations in DART, scenario **B** can be considered as most operational spruce forest representation, because the transformation of 3D shoots in turbid medium reduces considerably the simulation time and hardware requirements.

**Keywords:** terrestrial laser scanning; radiative transfer modelling; Norway spruce; 3D tree reconstruction; DART

<footnote n="1">¹ Global Change Research Institute, Czech Academy of Sciences, Bělidla 986/4a, 603 00 Brno, Czech Republic; homolova.l@czechglobe.cz (L.H.); hanus.j@czechglobe.cz (J.H.)</footnote>

<footnote n="2">² School of Technology, Environments and Design, Discipline of Geography and Spatial Sciences, University of Tasmania, Private Bag 78, Hobart 7001, Australia; zbynek.malenovsky@gmail.com</footnote>

<footnote n="3">³ CNRS, CNES, IRD, Center for the Study of the Biosphere from Space, Paul Sabatier University, 18 Avenue Edouard Belin, 31401 Toulouse, France; nicolas.lauret@cesbio.cnes.fr (N.L.); jean-philippe.gastellu-etchegorry@cesbio.cnes.fr (J.-P.G.-E.)</footnote>

<footnote>* Correspondence: janoutova.r@czechglobe.cz; Tel.: +420-511192-295</footnote>

<footnote>† Current address: Bělidla 986/4a, 603 00 Brno, Czech Republic.</footnote>
# 1. Introduction

Forests play a key role in the global environment as hotspots of terrestrial biodiversity [1,2], sinks of atmospheric carbon [3], and providers of multiple ecosystem services for human societies [4,5]. Green foliage is the main interface between a tree and its above-ground environment. It is involved in several key ecosystem processes, such as photosynthesis triggered by light interactions (i.e., absorption and scattering), as well as gas and water exchanges [6]. To understand the process-based functional and structural forest changes is a long-term scientific objective [7], and it can be supported by a multitude of optical remote sensing (RS) data and methods e.g., [8–14].

Optical RS methods require good understanding of the interactions between incoming solar radiation and the vegetation media that are dominated by photosynthetically active foliage. Canopy structure has a particular influence on canopy reflectance characteristics, and especially in the near and short-wave infrared and microwave optical domains [15–17]. Although light–canopy interactions can be studied effectively using physically based leaf and canopy radiative transfer models (RTMs), all RTMs rely on a certain level of simplification of vegetation canopy representations. Models simulating vegetation as a formation of horizontally homogeneous layers, such as the family of SAIL models by Verhoef [18] and Verhoef and Bach [19], the semi-discrete model of Gobron et al. [20], and the flux decomposition model of Kallel [21], are better suited to such structurally homogeneous canopies as those of crops, grasslands, or even tree crowns under the assumption of very high canopy closure and foliage density. Better suited to spatially heterogeneous and structurally complex forest canopies, however, are 3D RTMs working with tree architectures created from geometric shapes, such as the Flight model of North [22], FRT model of Kuusk and Nilson [23], or FLIES model of Kobayashi and Iwabuchi [24]. Steadily increasing computational capabilities allow for radiative transfer simulations of highly complex 3D canopies wherein the level of detail reaches to the individual leaf or needle. Such examples include the Raytran model of Govaerts and Verstraete [25], radiosity-graphics model of Qin and Gerstl [26], scalable plant RT model of Bailey et al. [27] and Discrete Anisotropic Radiative Transfer (DART) model of Gastellu-Etechgorry et al. [28,29].

In this study, we used the DART model, which simulates in-situ and RS acquisitions (LiDAR data and images in the optical and thermal spectral domains) and also the radiative budget of any landscape configuration, including terrain topography and atmosphere above the simulated scene [29]. Currently, DART enables constructing representations of trees by (i) creating simple predefined geometric shapes internally, (ii) importing external user-defined 3D objects, or (iii) transforming the imported 3D objects into simplified geometrical shapes [30]. The first of these utilities DART’s intrinsic functionalities defining a tree trunk and branches as simple 3D representations. Tree crown is represented by geometric shapes (e.g., an ellipsoid or a truncated cone) that are filled with triangles representing leaf facets or are created by cells known as voxels (i.e., 3D pixels) filled with a turbid medium representing foliage (Figure 1a). The turbid medium represents an infinite number of infinitely small leaves. The triangular facets and the turbid medium are characterized by predefined leaf optical and structural properties. Structural properties are primarily controlled by leaf area index (LAI) and leaf angle distribution (LAD), both specified either per tree or per canopy of the entire simulated scene. Other controlling parameters are distribution of empty voxels (i.e., air gaps) and leaf volume density inside a voxel, which can be specified either per crown or according to crown vertical level. Such a tree representation is computationally very effective and well suited for forest canopies with relatively simple tree architecture, such as those of broadleaf trees. This parametrization may be less suitable for conifers, however, due to its insufficient description of fine-scale scattering properties of needle shoots representing the main foliage elements. Therefore, we hypothesize that the second manner of constructing tree representations, in this case from user-defined 3D objects (Figure 1b), might be more appropriate for conifers and allow for producing more accurate DART simulations. The approach imports trees as groups of geometrically explicit 3D objects with different optical properties (i.e., trunks, branches, needle shoots of different age classes, etc.). DART can be executed directly over the original 3D objects constructed from triangular facets, which represents structurally the most accurate solution.
However, this option can be computationally highly demanding, especially when the imported 3D objects contains a lot of triangular facets. A third option is to transform the imported 3D foliage objects by voxelating them into leaf turbid cells (Figure 1c). Although such transformed crown representations are structurally less precise, they preserve all original aspects of crown shapes (e.g., distribution of air gaps) and foliage distributions (i.e., LAI and LAD). The voxelated nature of a forest stand representation substantially accelerates the computation and subsequently makes DART more operational with regards to its inverse modelling.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Three approaches for constructing trees using Discrete Anisotropic Radiative Transfer (DART): (a) creation of geometric crown shapes (basic representation) in DART, (b) import of external realistic 3D objects, and (c) transformation of realistic 3D tree objects into turbid-cells (voxelated) representations.</figcaption></figure>

Accurate physical representations that can account for complex 3D structures of ecosystems and link RS observations with studies of natural ecosystem processes are highly desirable. There is increasing demand, therefore, also for realistic 3D plant representations with natural architecture and foliage distribution. One group of methods for constructing 3D plants uses a set of parametric growing rules related to plant anatomy, topology, growth, and environmental factors derived from empirical observations. These are grammar-based methods, such as L-Systems [31], Grogra [32], OpenAlea [33], and recently developing computer graphic software (e.g., OnyxTree, Xfrog, and many others summarized at http://vterrain.org/Plants/plantsw.html). For instance, Widlowski et al. [34] used different virtual forests generated with Xforg and Arbaro software for the fourth phase of the radiative transfer model intercomparison (RAMI-IV) exercise. Nevertheless, most software tools create trees suitable only for exterior visualization and simulation of virtual landscapes (e.g., for computer games). They are less suitable for ecosystem studies requiring "biologically" correct 3D plant representations that account also for irregularities caused by natural and anthropogenic factors (e.g., spruce needle defoliation and regeneration; [35]).

Another group of methods for constructing 3D tree representations uses the rapidly developing technology of terrestrial laser scanning (TLS, also known as terrestrial LiDAR), which captures well the explicit 3D vegetation architecture [36–38]. Several studies have reported reconstruction of tree wooden skeleton (stem and branches) using TLS data by constructing a quantitative structure model (QSM) of topologically and hierarchically ordered components [39–41]. The TLS-based wooden skeleton reconstruction prefers data collected during leaf-off conditions, which means it is less than perfectly suited to conifers and evergreen-broadleaf mixed forests. In cases of evergreens, a TLS point cloud must be separated into wooden and foliage parts using (i) geometric properties [42], (ii) radiometric properties [43] or (iii) a combination of the two [44].
Considerably fewer studies have reconstructed a whole tree by adding foliage elements to a wooden skeleton. These studies concern mainly broadleaf tree species [45–47], and very few such attempts have been made for conifers [48]. Statistical information about spatial density, size distribution, and angular orientation of flat leaves could be extracted from the point cloud data [46,49,50]. For instance, Åkerblom et al. [46] introduced an algorithm to insert primarily broadleaf geometric 3D representations into QSM, but its extension to complex needle leaf objects considerably increased the computation requirements. Côté et al. [48] introduced the L-Architect model for conifers to reconstruct branching structure and to add foliage based on light availability derived from the point cloud data. Nonetheless, addressing 3D structure of conifers remains a challenge due to difficulties involving (i) leaf-wood separation, (ii) self-occlusion effects, and (iii) the complex geometrical nature of needle shoot elements.

The objectives of this study, therefore, were to (i) design an algorithm for reconstructing a close-to-reality 3D representation of Norway spruce trees from TLS measurements (i.e., reconstruction of wooden skeleton and natural distribution of needle foliage) and (ii) evaluate the effect of canopy structure on forest reflectance as simulated in a radiative transfer model. The second objective was achieved through DART simulations of four scenarios using different levels of spruce canopy structural details. The results of the scenarios were cross-compared and assessed against corresponding airborne hyperspectral and Sentinel-2 MSI RS reflectance observations.

## 2. Materials and Methods

### 2.1. Study Sites and Input Data

Computer reconstruction of trees, parameterization of DART, and the final assessment of the simulated canopy reflectance using airborne and satellite RS images required multiple data sources. All required data for this study were collected at two Norway spruce (Picea abies, L. Karst) forest sites located in the Czech Republic. The primary study site is located at the experimental research station Bílý Kříž in the Moravian–Silesian Beskydy Mountains (49°30′ N, 18°32′ E, 850–900 m a.s.l.). The site facilitates long-term research of tree ecophysiology and carbon fluxes. It is part of the national Czech Carbon Observation System (CzeCOS) and the international Integrated Carbon Observation System (ICOS) research networks e.g., [51,52]. The targeted stand is a spruce plantation growing on a southerly oriented slope (mean slope of 12.5°). Its characteristics in 2016 (when field and RS data were acquired) were as follow: mean tree age of 40 years, mean tree height of 17 m, mean stem DBH of 20 cm, mean stand density of 1230 trees ha⁻¹, and mean leaf area index (LAI) of 8.6 m² m⁻². The following additional stand characteristics were collected during a field and remote sensing campaign in August 2016 [53]: (i) foliar biochemistry and optical properties of needle samples for parameterization of DART simulations, and (ii) RS data (airborne hyperspectral and Sentinel-2 MSI images) used for evaluation of the DART reflectance simulations. Moreover, forest structural data on angular distribution of spruce shoots were obtained from an earlier study of Barták [54,55], who investigated 35-year-old spruce trees growing in the vicinity of the Bílý Kříž site.

Terrestrial laser scans (TLS) for 3D tree reconstruction could not be obtained for the Bílý Kříž site as the high density of forest stands did not allow collecting good TLS data on individual trees. Therefore, TLS data were acquired at a secondary study site, known as Černá hora and located in the Czech Republic’s Šumava National Park (48°58′22.4″ N, 13°32′53.8″ E, 1250–1300 m a.s.l.). The spruce stand at Černá hora (100-year-old, mean height of 25 m, stand density of about 800 trees ha⁻¹, SSE-oriented slope with gradient of 7°) had been severely damaged by multiple bark beetle outbreaks just 1–2 years before the TLS data acquisition. Therefore it was possible to acquire high-density laser scans of surviving individual trees that previously had been growing inside a canopy.
### 2.1.1. Forest Structure Data

Three-dimensional reconstruction of virtual spruce trees required two types of forest structural data: (i) TLS data for reconstruction of woody skeleton and spatial distribution of foliage (clumping) and (ii) data on geometrical and angular distribution of spruce needle shoots.

The TLS data were acquired using an OPTECH Ilris-36D (Teledyne Optech, Vaughan, ON, Canada) scanner at the Černá hora site in October 2009. The Ilris-36D TLS emits laser beams at the wavelength of 1.5 µm with a maximum field-of-view window of 40° × 40°. It records the amount and intensity of returned laser pulses at up to 2500 points per second. The laser scanning data were acquired for 15 mature spruce trees that had been left standing after clearing trees damaged by multiple bark beetle outbreaks. Each tree was scanned from two opposed positions, recording the first and the last echo return with the point sampling distance of 2 cm. The point clouds of the first and last return from the two scanning positions were aligned and co-registered in a single geometric coordinate system using the PolyWorks IMAalign module (InnovMetric Software Inc., Quebec City, QC, Canada). The PolyWorks IMSurvey module was used to filter out individual tree point clouds and to separate wooden structures from foliage points based on the intensity of the reflected photon signals. Considering that green foliage reflects about 10% and woody elements up to 50% of the laser signal emitted at 1.5 µm, we determined manually the intensity thresholds between foliage and wooden elements for each scan (first and last returns from two opposite scanning positions) separately. The intensity thresholds varied between 8% and 18% depending on a scan. The extracted point cloud of wooden elements was subsequently used to reconstruct a 3D tree skeleton consisting of trunk and main branches, whereas the foliage point cloud was used in the algorithm for biologically correct distribution of needle shoots within a tree crown (as described in Section 2.2).

It was too difficult to derive the angular distribution of coniferous shoots within a crown from the TLS data. Therefore, we extracted data from a destructive study of Barták [54,55], who had measured the inclination of spruce shoots in respect to their age and position within a crown (see Figure 2). Field sampling took place at the Bílý Kříž site in 1990–1991 and twenty 35-year-old spruce trees were analysed and we assumed that the measured inclination angles are representative for the current forest stand of the similar age.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Average angular distribution of spruce shoots according to their ages (c—current shoots, c-7 are 8-year-old shoots) and vertical position as measured by Barták [54] for 35-year-old Norway spruce trees at Bílý Kříž site.</figcaption></figure>
### 2.1.2. Needle Optical Properties

Optical properties (i.e., directional–hemispherical reflectance and transmittance) of spruce needles were measured during a field campaign at the Bílý Kříž site in August 2016. Needle samples of different age categories (current-, 1- and 2-year old needles) were collected from a sunlit and shaded branch of three trees growing at the main study plot in the vicinity of the flux tower. Needle optical properties were measured using a FieldSpec 4 spectroradiometer (Analytical Spectral Devices, Boulder, CO, USA) combined with an ASD RTS-3ZC integrating sphere using the method of gap fraction corrections suggested by Mesarch et al. [56] and further adjusted by Malenovský et al. [57] and Yanez-Rausell et al. [58]. Average needle spectra were used to parameterize the forest virtual scenes in DART, as described in Section 2.3.

### 2.1.3. Remote Sensing Data

Two types of RS data–airborne hyperspectral images and Sentinel-2 satellite data–were used for verification of DART-simulated canopy reflectance. Both RS images were acquired over the Bílý Kříž study site on 31 August 2016, close to the solar noon. Airborne hyperspectral images were acquired with a CASI sensor (Itres Research, Calgary, AB, Canada) in the range of 372–1044 nm (72 bands with spectral sampling distance of 10 nm and pixel size of 1.0 m) and with a SASI sensor (Itres, Research, Calgary, AB, Canada) in the range of 957–2442 nm (100 bands with spectral sampling distance of 15 nm and pixel size of 2.5 m). Corrections of the hyperspectral images were carried out in a data processing chain established at the CzechGlobe institute [59]. Radiometric corrections were performed using the factory calibration coefficients in the post-processing software developed by the sensors' producer. Geometric corrections (i.e., image orthorectification and geo-referencing) were performed using the GeoCorr software provided also by the sensors' manufacturer. Atmospheric and topographic corrections were made using ATCOR-4 software [60]. Images from CASI and SASI were combined into a single image hypercube (372–2442 nm) with final pixel size of 2.5 m.

The Sentinel-2 MSI image was corrected for atmospheric and topographic effects using the Sen2Corr toolbox provided by the European Space Agency. Spectral bands with nominal pixel size of 60 m and wide spectral bands were discarded, so that only the B2 to B7, B8a, B11, and B12 bands were retained in the final multispectral image cube.

Mean canopy hemispherical–conical reflectance function and its standard deviation were extracted from the CASI–SASI and Sentinel-2 MSI images for a region of interest (100 × 100 m) over the spruce forest study site in the vicinity of the flux tower, where field samples for needle optical properties and LAI measurements were obtained. The size of the region of interest was chosen to cover at minimum 4 × 4 Sentinel-2 pure forest pixels.

## 2.2. Reconstruction of 3D Spruce Tree Representation

The TLS point clouds supplied the explicit information about the position and orientation of dominant woody structures (trunk and branches) as well as about the foliage distribution and density of individual scanned trees. Woody and foliage elements were separated first using an intensity threshold that was determined for each tree scan manually based on visualization of laser return intensity histograms. The TLS data was itself insufficient for deriving the angular distribution of spruce needle shoots (considered here as the elementary units). Therefore, this information was taken from an external study on Norway spruce tree architecture [54,55] Figure 2.

Reconstruction of a 3D spruce representation from the TLS point cloud data was executed in the three stages documented in Figure 3:

1. Spatial separation and reconstruction of the wooden skeleton including trunk and main branches using an automated algorithm developed by Sloup [41] (Section 2.2.1, Figure 3b,c).

2. Translation and scaling of the foliage point cloud and the reconstructed wooden skeleton to fit site-specific tree dimensions. As the TLS data had been acquired at the Černá hora site for taller
trees, a scaling of dimensions was required to fit tree sizes at our Bílý Kříž primary research site (Section 2.2.2, Figure 3d).

3. Distribution of shoots within the tree crown by combining the information about foliage density from the TLS data with the ancillary information on shoot angular distribution [54,55]. (Section 2.2.3, Figure 3e).

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Reconstruction of a 3D spruce tree representation from terrestrial laser scanning data: (a) original terrestrial laser scanning (TLS) scan of a spruce tree, (b) point cloud separated into wooden elements, (c) reconstructed wooden skeleton, (d) reconstructed wooden skeleton with a scaled foliage point cloud, and (e) final 3D spruce tree representation populated by shoots of two age categories (current-year shoots in light green, earlier-year shoots in dark green).</figcaption></figure>

### 2.2.1. Reconstruction of Wooden Skeleton—Trunk and Main Branches

The algorithm for reconstruction of tree wooden skeleton had been developed for Norway spruce trees by Sloup [41] using TLS data acquired at the Černá hora site. It uses the method of Verroust and Lazarus [61] for extracting skeletal curves. It is fully automated and able to process point cloud data containing gaps from the omission of laser returns due to shadowing of any obstacles. It identifies spatially-related clusters of points (components) and the process of connecting branch components preserves the real architecture of branches and their connections to the trunk. The algorithm reconstructs the trunk and main branches, small branches at the crown periphery could not be reconstructed, because of insufficient details captured by TLS scans. Small wooden parts such as twigs can be included in the shoot representation.

The algorithm can be summarized into the following three steps (details and the code can be obtained from https://is.muni.cz/th/325196/fi_m/?lang=en):

1. Component identification: spatially related clusters of TLS points are identified and separated.

- constructing a neighbourhood graph to connect close points [61];
- the result is a set of components, each consisting of spatially related points belonging to the same branch or trunk part.

2. Component analysis: branch-like structures are reconstructed for each identified component.

- constructing a geodesic graph for each component that describes how the distance from the designated source point gradually increases throughout the component [61];

- collapsing points with similar distance from the source point, unless they belong to a different sub-branch.

3. Component connection: components are interconnected to form final trunk-branch structures.

- initial estimation of a trunk by fitting a segmented conical model to automatically selected trunk components;

- iterative selection of components to be connected to any other component while minimizing a cost function (combining the euclidean distance of the endpoints, outgoing vectors describing directions and probabilities of branch continuations, and slope function that roughly estimated expected slope of the main branches).

### 2.2.2. Translation and Scaling of Foliage Point Cloud

A reconstructed 3D wooden skeleton was built in a local coordinate system with its origin at the centre of the trunk base and neglecting original coordinate system of point cloud. A translation matrix between the two was established to ensure the same coordinate system for the foliage TLS point cloud as for the wooden skeleton. The reconstructed wooden skeleton and foliage point cloud were subsequently scaled to match site-specific tree dimensions. For this purpose, a scaling factor, $f_s$, was defined as the ratio between the required tree height and the maximum z-coordinate of the reconstructed wooden skeleton.

In our case study, the TLS data had been obtained from trees older than the spruces at our primary research site. Despite the difference in tree age of almost 60 years between the Černá hora and Bílý Kříž study sites, the discrepancy in the mean canopy height is not so very pronounced, as both forests are growing in mountainous environments and are exposed to similar environmental stress impacts. We therefore assume their woody and foliage architectures to be similar and possibly transferrable from one site to the another one. Scaled tree dimensions (crown diameter and height) were compared with existing forest allometry data showing good agreement. In addition, only foliage distribution is derived from the TLS point cloud data, not the actual density and angular distribution (described in the next section).

### 2.2.3. Shoot Distribution within a Crown

The third stage of tree reconstruction consists in distribution of the main foliage elements (i.e., needle shoots) within a tree crown. Because biochemical compounds together with the spatial and angular distribution of shoots influence light absorption and scattering, we considered in our reconstruction two needle age categories: (i) current-year needles, growing mainly at the crown periphery, and (ii) older needles, growing mainly inside the crown. The shoot distributing algorithm uses information about shoot positions from the foliage point cloud (a higher density of laser points indicates a higher density of leaves). The actual shoot density is controlled by the user through the definition of tree LAI. The leaf angle distribution (LAD) was adopted from the measurements of Barták [54] Figure 2. That algorithm has three main steps: (i) calculation of shoot positions, (ii) division of shoots into two age classes, and (iii) placement of shoots at predefined positions and angular orientations.
Calculation of Shoot Position

The number of shoots is given by the user-defined tree LAI. The total number of shoots within a reconstructed tree crown is computed as:

$$n_{shoots} = \frac{2 \cdot LAI \cdot S_{proj}}{S_{shoot}}, \qquad (1)$$

where $n_{shoots}$ is the total number of shoots within a crown, $S_{shoot}$ is total area of needles in one shoot (e.g., as calculated in the Blender software, www.blender.org), and $S_{proj}$ is the circular projection area of the tree crown represented by the extracted and scaled foliage point cloud. The diameter $d$ of the circular projection area $S_{proj}$ is calculated as:

$$d = \min(\max_{x} - \min_{x}, \max_{y} - \min_{y}), \qquad (2)$$

where $\max_{x}, \max_{y}, \min_{x}$, and $\min_{y}$ are the maximum and minimum coordinates of the foliage points in the $x$ and $y$ dimensions. Calculation of the shoot positions uses a grid of voxels. The foliage point cloud is, therefore, first separated into a grid of voxels of a given size (in our case 0.18 m). The voxel size influences the computational time, as placing more points within a single voxel requires more time for computing all shoot positions. Hence, the voxel size should be user-specified according to the computational resources available and the density of the foliage point cloud.

The actual shoot positions are calculated using the $k$-means function based on the cluster analysis theory. Each voxel contains a set of points $\mathbf{x} = \{\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_n\}$, where $n \in \mathbb{N}$ is the number of points inside the voxel currently being processed. Each point is defined by its position $\mathbf{x}_i = (x, y, z)$. The $k$-mean clustering partitions points into $k < n$ subsets $\mathbf{S} = \{S_1, S_2, \dots, S_k\}$. The $k$ number of subsets is then computed from the total number of shoots within the crown ($n_{shoots}$, Equation (1)) as:

$$k = \frac{n}{\frac{n_{points}}{n_{shoots}}}, \qquad (3)$$

where $n_{points}$ is the total number of points in the whole foliage point cloud. The $k$-mean function computes also the coordinates of the cluster centroid, which then represents the actual shoot position. The algorithm iterates through all voxels containing at least one laser point.

Separation of Shoots by Age Categories

For this step, we need to identify a tree crown envelope to assign current shoots at the periphery and older shoots inside the crown. Therefore the foliage cloud point has to be divided into a voxel grid again and we recommend setting the size of voxels at approximately twice larger than that for the shoot positioning grid (in our case we use the voxel size of 0.3 m). The tree crown envelope and its width were defined using an iterative approach that is presented in Appendix A. The tree was then divided into several vertical layers ($l_n$) and a percentage of current-year shoots was assigned to each layer based on the field measurements, as summarized in Table 1 (details in Appendix B). Different optical properties (Figure 4) were defined for each needle category in order to capture the natural biochemical variability.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Measured optical properties of Norway spruce needles and other stand elements used for parameterization of DART simulations: (a) optical properties for the top, middle, and bottom crown levels applied in the case of basic tree representation; (b) optical properties of two needle age categories (current-year and older needles) that were used for the 3D tree representation reconstructed from TLS point clouds; and (c) optical properties of spruce bark and forest floor.</figcaption></figure>

<table><caption>**Table 1.** Field data on spatial and angular distribution of Norway spruce shoots (destructive sampling of 20 trees) extracted from Barták [54,55]. Percentage of current-year shoots, mean shoot elevation angle, and that angle's variation in the vertical crown levels associated with branch whorls are counted from top to base.</caption><thead><tr><th>Whorl Order</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>&gt;10</th></tr></thead><tbody><tr><td>whorl height (m)</td><td>0–1</td><td>1–2</td><td>2–3</td><td>3–4</td><td>4–5</td><td>5–6</td><td>6–7</td><td>7–8</td><td>8–9</td><td>&gt;9</td></tr><tr><td>% of current shoots</td><td>1</td><td>0.82</td><td>0.61</td><td>0.49</td><td>0.35</td><td>0.29</td><td>0.23</td><td>0.17</td><td>0.07</td><td>0.05</td></tr><tr><td>mean shoot elevation angle (°)</td><td>35</td><td>35</td><td>5</td><td>–15</td><td>–35</td><td>–45</td><td>–25</td><td>–35</td><td>–32</td><td>–35</td></tr><tr><td>elevation angle variation (°)</td><td>30–40</td><td>30–40</td><td>0–20</td><td>–20–15</td><td>–45–5</td><td>–50––5</td><td>–40––5</td><td>–50––20</td><td>–50––20</td><td>–50––20</td></tr></tbody></table>

## Shoot Types and Their Positioning within the Crown

Because needle shoots can be constructed as 3D objects, their shape and complexity can be altered according to user and computational requirements. Here, we used two different 3D shoot models: (i) a shoot constructed from individual non-flat needles along a central twig (Figure 5a), and (ii) a shoot simplified to radial planes with total area equal to the leaf area of single needles while maintaining the original needle lengths (Figure 5b). The main difference between the shoot models consists in the number of object facets. The shoot model with individual needles contained 824 facets per shoot, whereas the simplified model with planes contained only 16 object facets. Both shoot representations had been created and used for the RAMI-IV radiation transfer model intercomparison exercise [34].
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Example of a 3D representation of a spruce shoot with individual needles (a), and a simplified representation with planes along the central shoot axis (b) having total area equal to the leaf area of the corresponding needles.</figcaption></figure>

Once the exact positions for shoots have been determined, shoots as individual 3D objects are distributed within the tree crown. The shoots must first be rotated according to predefined angular boundaries varying in the vertical direction as summarized in Table 1 and Figure 2. Finally, shoots are translated to their actual positions inside the crown. A transformation matrix for angular rotation and spatial translation is calculated for each shoot separately, as described in Appendix C.

## 2.3. DART Forest Scene Parametrization

Small forest scenes of 10 × 10 m and with a fixed distribution of 10 reconstructed trees were created in DART using either internal basic tree shapes filled with a turbid medium (Figure 1a, Appendix D) or the external 3D representations reconstructed from TLS point clouds (Figure 1b,c). Main differences among the parameterization approaches are not only in the tree structure descriptions but also in definitions of the leaf optical properties, LAD, and LAI. The remaining DART input parameters were tailored to forest conditions captured at the Bílý Kříž study site during the field and remote sensing campaign in 2016. These include tree dimensions; optical properties of forest floor, spruce bark and needles; actual sun position during the acquisition of airborne and Sentinel-2 images; and description of the atmosphere (based on a mid-latitude summer model in Modtran [62]). The input parameters common to all DART simulations are summarized in Table 2, and the optical properties are depicted in Figure 4c.

<table><caption>**Table 2.** Basic DART scene parameters.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>scene dimension</td><td>10 m</td></tr><tr><td>cell dimension</td><td>0.5 m</td></tr><tr><td>terrain</td><td>flat</td></tr><tr><td>number of trees</td><td>10</td></tr><tr><td>scene canopy cover</td><td>80%</td></tr><tr><td>scene LAI</td><td>8.5</td></tr><tr><td>average tree height</td><td>15 m</td></tr><tr><td>sun zenith angle</td><td>44°</td></tr><tr><td>sun azimuth angle</td><td>153°</td></tr></tbody></table>

Leaf optical properties were computed differently for the trees created internally and externally. In the case of basic DART tree shapes, they were mixed based on the ratio of current-year and older shoots (Figure 4a), as follows: (i) 7:3 for a sunlit vertical crown layer (top 17% of tree crown), (ii) 1:1 for a transitional sunlit and shaded crown layer (middle 30% of tree crown), and (iii) 1:8 for a shaded crown layer (bottom 53% of tree crown). In the case of the reconstructed 3D trees, the leaf optical
properties were assigned according to the distribution of two needle age categories as described in Section 2.2.3 (Figure 4b).

## 2.4. Evaluation of the Spruce Forest RT Simulations

All RTM simulations were conducted with DART model version 5.7.1 (build 1061) according to four scenarios representing various levels of canopy structural complexity as summarized in Table 3. Scenario A, basic turbid, which is regarded as structurally the simplest, uses basic predefined crown shapes and represents the foliage as turbid cells. The remaining three scenarios, B–D, use 3D trees reconstructed from TLS, but they differ in the complexity of 3D shoot objects and in the foliage representation as object facets or as turbid cells transformed from these objects. Scenarios A–C were simulated for the spectral range from 450 to 2300 nm with the 2 nm bandwidth and subsequently convolved to the spectral resolutions of CASI–SASI and Sentinel-2 MSI images. Scenario D, needle facet, is structurally the most complex and thus computationally the most demanding scenario. Therefore, only 13 spectral bands with the central wavelengths of Sentinel-2 MSI bands were computed for this scenario.

<table><caption>**Table 3.** Four scenarios representing varying structural complexity of a virtual Norway spruce forest scene used to simulate canopy reflectance in the DART model.</caption><thead><tr><th>Scenario Label</th><th>3D Tree Rep.</th><th>3D Shoot Rep.</th><th>DART Forest Scene Types</th></tr></thead><tbody><tr><td>A basic turbid</td><td>DART predefined tree crown shapes</td><td>not applicable</td><td>turbid cells</td></tr><tr><td>B needle turbid</td><td>reconstructed from laser scans</td><td>exact 3D needles</td><td>transformed into turbid cells</td></tr><tr><td>C plane facet</td><td>reconstructed from laser scans</td><td>simplified 3D shoots as planes</td><td>3D shoots with object facets</td></tr><tr><td>D needle facet</td><td>reconstructed from laser scans</td><td>exact 3D needles</td><td>3D needle object facets</td></tr></tbody></table>

From DART-simulated images of each scenario we extracted mean scene spectral reflectance and its standard deviation, as well as mean reflectance of sunlit crown parts and its standard deviation. DART-simulated spectral signatures were compared with the canopy reflectance signatures extracted from CASI–SASI and Sentinel-2 images. Root mean square error (RMSE) and index of agreement *d* between simulated ($R_{SIM}$) and observed ($R_{OBS}$) spectral signatures were computed as:

$$RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (R_{i,SIM} - R_{i,OBS})^2}, \qquad (4)$$

$$d = 1 - \left( \frac{\sum_{i=1}^{n} (R_{i,SIM} - R_{i,OBS})^2}{\sum_{i=1}^{n} (|R_{i,SIM} - \overline{R_{OBS}}| + |R_{i,OBS} - \overline{R_{OBS}}|)^2} \right), \qquad (5)$$

where *n* is total number of either CASI–SASI or Sentinel-2 MSI spectral bands, $\overline{R_{OBS}}$ is mean value of observed spectral signatures. The CASI–SASI spectral bands that are strongly influenced by atmospheric water vapour absorption (i.e., 893–1002 nm, 1107–1167 nm, 1302–1527 nm, and 1737–2052 nm) were removed from the RMSE and *d* computation. Further, we extracted and cross-compared multi-directional reflectance distributions of the simulated canopy scenarios.

The main criteria to evaluate the optimum scenario are (i) overall agreement with real observations extracted from airborne hyperspectral and Sentinel-2 image data and (ii) computation requirements.

# 3. Results

## 3.1. Reconstructed 3D Spruce Trees

The total number of leaf facets of a single tree constructed according to TLS scans (Figure 6) varied between 11 million and 30 million, depending on crown dimensions (see Table 4). The LAI of each individual tree crown, i.e., relative to the tree crown vertical projection onto the ground, was set to 12 m² m⁻² (Equation (1)). This value had been found to be the optimal LAI for building up a forest
scene with canopy cover of approximately 80% and canopy LAI equal to 8.5 m² m⁻² (as measured in the field). The DART forest scene of 10 trees was built from the four individual tree representations, which were replicated, rotated, and re-distributed to achieve the desired canopy cover while avoiding a regular spacing (Figure 7).

<table><caption>**Table 4.** Structural characteristics of the four reconstructed 3D spruce trees (Figure 6) populated with complex needle shoot objects (leaf area index = 12 m² m⁻²).</caption><thead><tr><th>Parameter</th><th>Tree 1</th><th>Tree 2</th><th>Tree 3</th><th>Tree 4</th></tr></thead><tbody><tr><td>tree height (m)</td><td>15.0</td><td>15.1</td><td>16.3</td><td>14.0</td></tr><tr><td>crown length (m)</td><td>13.5</td><td>10.5</td><td>13.0</td><td>11.5</td></tr><tr><td>crown diameter (m)</td><td>2.3</td><td>2.5</td><td>3.5</td><td>3.7</td></tr><tr><td>number of wood facets</td><td>112, 448</td><td>106, 112</td><td>186, 432</td><td>197, 760</td></tr><tr><td>number of leaf facets</td><td>11, 617, 576</td><td>16, 829, 376</td><td>30, 385, 824</td><td>30, 791, 232</td></tr><tr><td>proportion of current shoot facets</td><td>0.40</td><td>0.47</td><td>0.29</td><td>0.22</td></tr></tbody></table>

Rigorous validation of the reconstructed 3D trees would be difficult and tedious, as it would require that the evaluated structural parameters of the reference trees be measured independently from their TLS scans. Moreover, the original TLS data were spatially scaled to different dimensions and the density of needle shoots was altered according to the user-defined value of tree crown LAI. Therefore, the distribution of shoots and the overall habitus of reconstructed trees were inspected only through visual comparison with the original TLS point clouds (see Figure 3a) and photographs of modelled trees. Later, we evaluated the reconstruction of 3D trees in the context of the forest stand canopy reflectance simulated in DART and compared to the real airborne and space-borne observations.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Four individual Norway spruce tree 3D representations reconstructed from TLS data point clouds using complex needle shoots. Colour rendering was done in the Blender software (www.blender.org), showing wood components in brown, newest needle shoots in light green, and older shoots in dark green colours.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Visualization of Norway spruce forest 3D scenes (a–c) and DART-simulated top-of-canopy reflectance images for scenarios A–D. (a) scenario A, *basic turbid*, with DART predefined crown shapes and foliage as turbid cells (solid green voxels represent foliage with predefined optical properties, whereas semi-transparent light green voxels are foliage intersected by 3D woody objects), (b) scenario B, *needle turbid*, with 3D trees reconstructed from TLS data occupied by turbid cells converted from facet-based needle shoots (solid brown voxels represent foliage with mixed optical properties of current-year and older-needle shoots), and (c) scenario C, *plane facet*, and scenario D, *needle facet*, with 3D trees reconstructed from TLS data having shoots approximated as planes or retained as needle facets. The upper images (A–D) show nadir false colour composites (red = 779 nm, green = 666 nm, and blue = 552 nm) of simulated forest canopy reflectance for each scenario.</figcaption></figure>

## 3.2. Simulations of Forest Canopy Reflectance in DART

DART simulations of Norway spruce forest canopy reflectance between 450 and 2300 nm (i.e., 926 bands with step and width of 2 nm) for scenarios A–C were carried out on a local server (40 cores, 350 GB RAM). Simulations of the turbid cell scenarios A and B were computed relatively quickly (each scene in less than 3 h, using only 2 processor cores), but to retain the foliage elements as object facets required substantially increased computation time. The forest scene of scenario C (shoots as planes) contained 3.7 million leaf facets and its computation took almost 34 h (using 8 processor cores). Computation of scenario D, *needle shoots*, containing 189 million facets, was extremely demanding computationally and had to be processed on the Czech national computer grid MetaCentrum (https://metavo.metacentrum.cz/en/index.html).

Figure 7 displays the reconstructed virtual forest scenes and their corresponding false-colour composite images of canopy reflectance simulated in DART. Scenario A, *basic turbid*, shows regular shapes and homogeneous filling of spruce crowns, which are not often seen on real airborne spectral images. The structurally more heterogeneous scenarios B–D revealed larger spatial variability of reflectance within and outside of the crown pixels. One can observe more shaded pixels inside and at the crown edges and also more sunlit forest floor pixels (the cyan colour), especially in scenarios C and D. Despite that crown diameters were the same in all scenarios, retaining the real branch architecture resulted in tree crown peripheries being more transparent when compared to scenario A.

Overall, the greatest simulated canopy reflectance in green and near-infrared (NIR) spectral domains was found in scenario A. Simulated red and shortwave infrared (SWIR) reflectance values were converging in all simulated scenarios, which indicates that these spectral wavelengths are less affected by differences in simulated canopy architectures (Figures 8). All simulations were carried
out with added sensor viewing directions around the hotspot (i.e., hotspot oversampling). One can observe that turbid media (scenarios A and B) produced a sharper hotspot feature compared to facet-based simulations (scenarios C and D) (Figure 8), which is an artifact of the turbid canopy hotspot analytical computation. Inspecting closely the multi-angular reflectance behavior of the four scenarios (Figure 9), one can note overestimations in the green wavelength (552 nm) for the turbid scenarios (A and B) when compared to the facet-based scenarios (C and D), especially at oblique zenith viewing angles. All multi-angular behaviors agree well for the red reflectance band (666 nm), as they were influenced by strong chlorophyll absorption. The only exception relates to the *plane facet* of scenario C, which indicates a reflectance overestimation in oblique back-scattering directions. The largest multi-angular reflectance differences can be seen in the NIR band (776 nm). The *basic turbid* scenario A overestimates, while the *plane facet* scenario C underestimates, NIR reflectance in all viewing directions when compared to scenario D. Finally, all simulated scenarios produced a consistent reflectance with minimal differences in the SWIR wavelength (1647 nm). Minor overestimations were observed for the *needle turbid* scenario B in strongly oblique (> 60°) forward-scattering viewing directions.

The first derivative (FD) of simulated canopy reflectance reveals changes in the shape of spectral signatures. FD was not calculated for scenario D, for which shape changes would be too abrupt, due to its inclusion of just 13 simulated spectral bands. Results for scenarios A–C show that, despite their differences in absolute reflectance, the signature shapes are very similar, especially for the NIR and SWIR wavelengths that are not impacted by atmospheric water vapour absorption (Figure 10). The largest differences were found in the green and mainly in the red-edge region. One can observe substantial FD peaks, which correspond to inflection points on the steepest parts of the reflectance curves. In most bands, the FD of scenarios A, *basic turbid*, and *B*, *needle turbid*, are similar, which likely is caused by their turbid cell nature. Interestingly, the reflectance signature shape closest to that of the real airborne data acquired over the study site using CASI and SASI sensors was found in the case of scenario C, *plane facet*.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** DART-simulated forest canopy reflectance in the solar principal plane for the four simulated scenarios (A–D) in four spectral bands (a) green—552 nm, (b) red—666 nm, (c) NIR—800 nm, and (d) SWIR—1647 nm.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Multi-angular distribution of simulated forest canopy reflectance in four spectral bands (green—552 nm, red—666 nm, NIR—779 nm, and SWIR—1647 nm) for the most complex canopy scenario **D**, *needle facet*, (left plots), followed by the reflectance differences between scenario **D** and other three scenarios (**A**, **B**, and **C**).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** The first derivative computed for nadir spruce canopy reflectance simulated for three scenarios **A**, *basic turbid*, **B**, *needle turbid*, and **C**, *plane facet*, by DART and extracted from the airborne imaging spectroscopy data acquired with CASI and SASI sensors. The missing parts of spectral signatures are those wavelengths strongly impacted by atmospheric water vapour absorption.</figcaption></figure>

### 3.3. Comparison with Airborne and Satellite Data

Simulated canopy reflectance signatures for the four scenarios (**A**–**D**) were plotted and compared against spectral signatures extracted from the corresponding airborne hyperspectral CASI-SASI images and satellite Sentinel-2 MSI image (Figures 11 and 12). Taking advantage of the high ground spatial detail of airborne hyperspectral images, we compared the simulated scenarios with (i) mean forest stand spectra extracted from CASI-SASI image with pixel size of 2.5 m and entire DART scenes
(Figure 11a,b), (ii) spectral signatures extracted from sunlit tree crowns of CASI and DART images of 1.0 m pixel size (without contributions from canopy shaded parts and background) (Figure 11c,d), and (iii) mean canopy spectral signatures extracted from Sentinel-2 MSI image with pixel size of 20 m and DART scenes (Figure 11e,f). The comparison of forest stand reflectance extracted from airborne CASI-SASI and Sentinel-2 MSI images for the same region of interest showed good agreement, even though the reflectance variability in the airborne data is naturally greater due to its finer pixel size (2.5 m) (cf., Figure 11a,e). The consistency between the two independent RS observations gives us a certain level of confidence in using the RS data as a robust benchmark for evaluating the different scenarios of spruce forest parameterization in DART. When compared to the RS observations, all scenarios (A–D) overestimated the stand canopy reflectance across the entire spectral region between 450 and 2300 nm. Scenario A, *basic turbid*, yielded the largest differences from the RS observations, mainly in the visible and NIR regions. Other implementations of the tree 3D architectures (scenarios B–D) brought the DART forest reflectance simulations closer to the RS observations.

In the visible region, all scenarios systematically yielded greater reflectance than did RS observations. The simulations manifest a systematic reflectance overestimation by 0.01–0.02 in green peak. The red reflectance of the chlorophyll-a absorption feature around 675 nm also exceeded that of the RS observation by 0.01, and in the case of scenario C, *plane facet*, this was by even 0.02. This deviation disappeared, however, when inspecting the reflectance of only sunlit parts of tree crowns (cf., Figure 11b,d). This suggests that the mean reflectance of scenario C might be influenced by sunlit background pixels that are substantially increasing simulated values in red wavelengths.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Comparison of DART-simulated nadir reflectance (BRF) signatures with (a,b) mean forest stand reflectance extracted from airborne hyperspectral CASI-SASI images (pixel size 2.5 m), (c,d) mean reflectance of sunlit tree crowns extracted from DART-simulated images and airborne hyperspectral CASI image (pixel size 1.0 m), and (e-f) mean canopy reflectance extracted from satellite Sentinel-2 MSI image (pixel size of 20 m). The left plots show close-up details of the visible spectral range (450–700 nm).</figcaption></figure>

The largest reflectance variability among the simulated scenarios and differences from RS data were obtained for the NIR spectral region. When compared to the airborne data, the NIR reflectance around 850 nm of scenario A, *simple turbid*, was as much as 0.1 larger. This difference increased to 0.15 for sunlit crown spectra. A similar NIR reflectance overestimation was observed also for scenario D, *needle facet*. NIR reflectance of scenarios B, *needle turbid*, and C, *plane facet*, fell within the standard deviation interval of the RS observations. Finally, all scenarios simulated very similar reflectance for the SWIR region. Nevertheless, all scenarios systematically overestimated the RS observations around 1650 nm by 0.04 and around 2200 nm by 0.02.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Scatterplots between DART-simulated and remote sensing (RS) data extracted from spruce forest stand reflectance. Top row: average stand reflectance extracted from CASI-SASI image, middle row: average reflectance of sunlit tree crowns extracted from CASI image, and bottom row: canopy reflectance extracted from Sentinel-2 MSI image (only narrow spectral bands with spatial resolution ≤20 m that are not affected by atmospheric water absorption are plotted). In order to provide statistically comparable data sets, the open circles in the case of airborne data (top and middle row) represent wavelengths corresponding to the Sentinel-2 MSI spectral bands.</figcaption></figure>

Scatterplots, similarity indicator *d*, and *RMSE* between the DART-simulated scenarios and RS observations (airborne and Sentinel-2 MSI) were plotted and computed only for wavelengths corresponding to those Sentinel-2 MSI spectral bands not affected by atmospheric water absorption, (i.e., not interpolated during the atmospheric corrections) (Figure 12). The highest *d* and lowest *RMSE* values were obtained for scenario C, *plane facet*, followed by *B*, *needle turbid*, and *D*, *needle facet*.

# 4. Discussion

## 4.1. Reconstructed 3D Spruce Trees

The first goal of this study was to create a biologically correct 3D digital representation of Norway spruce trees for purposes of forest canopy radiative transfer modelling [63,64]. The results demonstrate that the designed algorithm is capable of reconstructing the architecture and geometrical properties of spruce wooden and crown foliage elements from TLS observations of individual trees. Quantitative validation of the reconstructed 3D spruce representations could not be performed, as doing so would require reference measurements of tree structural (geometrical) features and the scanning device position, which were not at our disposal. Such a detailed quality assessment is not even needed, however, because the aim was not to reconstruct a given tree as accurately as possible but rather to obtain biologically correct 3D tree objects that can form a representative forest canopy with a variable leaf area index and user-defined foliage biochemical traits. Consequently, the developed method allows one for reconstructing trees with any required value of LAI in order to match the forest stand LAI measured at the site where the RS data was acquired.
The algorithm is automated and works with separated point clouds of wooden skeleton and foliage elements. The fidelity of a reconstructed tree architecture depends on the quality of TLS scans, which drives the decision on either manual or automatic separation of wood and foliage point clouds. The wooden skeleton reconstruction algorithm [41] correctly handles point gaps in wood point cloud data caused by occlusion of emitted laser beams due to branches within the scanner field of view. The algorithm for foliage distribution has been designed specifically for Norway spruce trees, but it might be applicable also to other spruce species and, after necessary adjustments, to other tree biological species. A user has the possibility to control tree LAI by adjusting the density of foliage objects (in this case shoots) in proportion to the density of TLS points. The angular orientation of spruce shoots had to be set according to manual field measurements of a limited sample of trees [54], as it was not feasible to retrieve this information from TLS scans. The leaf angle distribution might be obtained for broadleaf trees e.g., [49], but the method is not applicable for conifers with small narrow needles organized in shoots.

The approach developed in this study can realistically reconstruct the architecture of an entire spruce tree, including wooden skeleton and foliage geometry, while also implementing different groups of needle shoots according to their age. This ability is crucial for such coniferous tree species as Norway spruce, where leaves or shoots emerging in different vegetation seasons (i.e., several generations) co-exist within a single tree crown. Previous studies (e.g., [63–66]) solved this problem by dividing tree crowns into several vertical layers, to which different needle optical properties are attributed according to measured needle-age ratios and sun radiation availability. Our approach allows one for populating a tree skeleton with real needle shoots, thus reflecting the realistic distribution of specific leaf optical properties according to their age. This capability is especially important for simulations of high spatial resolution optical images (with pixel-size on the order of centimetres), where diversity in both, optical properties and forest structure has an important influence on the forest canopy reflectance. In this respect, the reconstruction of wooden skeleton and crown architecture might be sufficient, but the distribution of optical properties could be treated more realistically. Besides shoots, young branches and younger parts of the trunk are in the case of spruce also populated by needles. In our study, however, all branches and trunk parts were assigned with the optical properties of bark. Implementing needles, or at least their optical properties, for the top parts of branches and trunks might mean a substantial difference in simulated red and NIR canopy reflectance, especially for the nadir viewing direction.

Another possible improvement is to create a more realistic reference needle shoot 3D model. The representation used in this study was relatively simple. Its design was ruled by the limited capabilities of handling 3D objects in older versions of DART and limited availability of computation resources. Moreover, a universal distribution of needles along the central wooden twig does not resemble the actual variability existing in natural shoot geometrical formations. The architecture of each shoot varies depending on illumination conditions, age, and stress impacts. More genuine 3D shoot representations could be obtained from high-resolution laboratory laser scanning. Implementing several types of shoot representations reflecting the natural variability driven by shoot age and location within the tree would improve the accuracy of the radiative transfer modelling, and particularly that of the light scattering within shoots.

## 4.2. *Simulations of Forest Canopy Reflectance in DART*

One of the key issues for remote sensing of coniferous forests is the upscaling of optical properties from needles to canopies. From an operational point of view, it is easier to measure optical properties of needles rather than of shoots. Inasmuch as shoots are the main scattering elements in coniferous canopies, however, there is a strong demand for characterizing the shoot optical properties and phase functions with high accuracy [67–70]. Traditional forest reflectance models based on the radiative transfer equations handle the shoot-level needle clumping by correcting the radiation attenuation coefficient with a clumping index. That compensates for a reduction in the interception of radiation by
the canopy at a fixed LAI [68]. Hence, implementation of the precise 3D spruce representations can help in disentangling and better understanding light interaction related to the needle foliage clumping at different structural levels (i.e., shoots/twigs, branches, crowns, and stands).

The reflectance simulation of whole forest canopy in DART served to compare different approaches of tree parametrizations and to assess their performance in terms of accuracy, operational implementation, and feasibility. To produce DART simulation of a forest canopy scene 10 × 10 m in size and with 3D tree objects created from almost 189 million facets was until recently unfeasible. We were able to execute such a geometrically detailed simulation using high-performance computation resources of the Czech Republic's national grid organization known as MetaCentrum. Even then, we were able to simulate only 13 spectral bands and total computation time took more than three weeks. Outputs of this most-detailed simulation were used in comparison with other, simpler parametrizations, where we are able to simulate hyperspectral images with the 2 nm spectral resolution in much shorter time and using substantially fewer computational resources. With future DART code improvements and availability of more powerful computation resources, we expect to simulate full hyperspectral images even for the most detailed 3D representations, as in the scenario **D**.

The results confirmed an expectation that simulations using the reconstructed 3D representations produce more consistent canopy reflectance, especially in red and NIR wavelengths (see Figure 8). The first derivative of the canopy reflectance revealed that the most substantial shape changes in the reflectance function are between the *basic turbid* scenario A and the other scenarios at the red edge spectral region (Figure 10). Nevertheless, depending in each case upon the purpose of their use, the results of all scenarios could be found acceptable. Scenario A, *basic turbid*, differs most from the others. It is the least demanding of computational resources and its reflectance results are in agreement with the reference scenario **D** in the red and SWIR wavelengths. Contrary to this scenario **D**, *needle facet* requires high-performance computational resources and a substantial simulation time, due to detailed 3D needle shoots. Reflectance results of scenario **B** *needle turbid* are in most cases very similar to the most detailed scenario **D** *needle facet*. The only substantial difference was found in the NIR region (800–1200 nm). As expected, canopy reflectance of both scenarios simulating canopy as a turbid medium behaved differently in the hot spot region, due to the turbid simulation's hot-spot analytical solution. Finally, compared to scenario **D**, the scenario **C** *plane facet* slightly underestimated red reflectance, especially around the hot-spot region, and overestimated NIR reflectance in all simulated viewing directions (see Figure 9).

### 4.3. Comparison with Airborne and Satellite Data

Widlowski et al. [34] pointed out that the reconstruction of a natural forest stand scene for purposes of modelling radiative transfer could be problematic, because several tree habitus characteristics might be unavailable or difficult to measure. Many canopy biochemical, optical, and mainly structural properties, such as the geometrical distribution of foliage elements (leaves or shoots), in a canopy vertical profile or allometric parameters of wooden elements, including the branching angles, are often not measured simultaneously with RS spectral image acquisitions. The Bílý Kříž research station in Czech Republic is an ecological monitoring site where ground measurements of the monitored spruce stand were collected simultaneously with airborne imaging spectroscopy data (CASI and SASI sensors) and satellite Sentinel-2 MSI images. This allowed us not only precisely to parameterize the DART 3D scene, but also to compare the simulations with RS data acquired during the same vegetation phenophase and at the same location.

The comparison revealed a systematic overestimation of simulated bidirectional reflectance factors, especially within the spectral interval of 450–570 nm. This mismatch could be caused by a greater contribution of background, represented in our simulations by needle litter having relatively high visible spectral region of reflectance. The possibility that this is the cause is supported by results obtained from separating spectral responses of sunlit crowns within the simulated and acquired images. Suppression of the understorey influence in red wavelengths is obvious in the case of scenario
C, *plane facet*, (cf. Figure 11b,d). Another possible cause of this discrepancy might be a mismatch in the distribution and content of shadowed pixels (between tree crowns) in RS data and simulations. The physical modelling treats light transfer according to optical laws without any noise, while the RS sensors and processing of their data introduce several types of noise that DART does not simulate. Although we used the needle optical properties measured during the RS data acquisitions, a certain level of non-representativeness originating from a disproportionate ratio between current and older needle generations might still occur and contribute also to this mismatch. Finally, a spectral shift between simulated and RS data around 1800 nm is caused by the interpolation of airborne data reflectance in the course of the atmospheric corrections.

The comparison of canopy reflectance extracted from RS data and DART simulations, as if mimicking narrow spectral bands of Sentinel-2 MSI with the spatial resolution ≤20 m and unaffected by atmospheric water vapour absorption (Figure 12), revealed almost the same results for the last three scenarios (B, C, and D). This suggests that the most efficient and optimal approach for simulating the Sentinel-2 MSI forest stand reflectance in DART is to use scenario B, *needle turbid*, which can deliver comparable reflectance simulations as scenarios C and D, but the computation time is substantially shorter.

**5. Conclusions**

In this study, we evaluated the influence of Norway spruce forest structural representations on DART modelled top-of-canopy reflectance. DART simulated nadir reflectance values were compared with actual airborne hyperspectral and satellite multispectral images acquired over the simulated research site. Because the multi-angular directional reflectance measurements of the investigated spruce forest were not available, the multi-angular simulations were compared against the most detailed and presumably most accurate scenario D, *needle facet*. The results revealed improvements in reflectance accuracy simulated for scenarios wherein the 3D forest scene representations were retrieved and parametrized using the method designed in this study. We identified scenario B, *needle turbid*, as the most effective scenario for nadir simulations. This scenario was able to produce the visible and NIR reflectance with an acceptable level of accuracy (i.e., similar to that of scenario D and close to actual observations) with relatively short computational time and without a need for high-performance computer resources. Being based on the *needle turbid* representation of the stand, however, scenario B could not very well reproduce reflectance around the hot-spot and in the extreme zenith angles. In the case of the multi-angular directional reflectance, scenario C, *plane facet*, tracks closely the most detailed scenario D, *needle facet*, while maintaining a comparable (in NIR) or even better match (in visible wavelengths) with the reflectance of RS data acquired from the nadir.

The 3D representation of Norway spruce trees reconstructed from TLS data with the semi-automatic algorithm designed in this study was found to reassemble adequately the habitus of real Norway spruce trees and increase the fidelity of canopy reflectance simulations in DART. The topographically accurate capturing of the tree structures and more detailed parametrization of foliage optical properties throughout a tree vertical profile allowed for more accurate solar radiation scattering within a simulated forest stand. Because the results obtained showed reflectance overestimations in visible and NIR spectral regions when compared with RS observations, there clearly is room for further improvements. Subsequent research should be focused upon (i) more realistic distribution of optical properties for wooden parts of forest stands (especially for young branches and tips of tree crowns), and (ii) the extension and improvement of needle shoot representations. Several needle shoot types could be reconstructed and used instead of a single, universal 3D model. Their size and geometrical parameters should reflect natural variability driven by shoot location within the tree crown, age, needle angular distribution, and within-shoot clumping in order to produce more realistic spruce tree representations for purposes of estimating radiative transfer.

**Author Contributions:** Conceptualization, R.J., L.H. and Z.M.; Data curation, R.J., L.H. and J.H.; Formal analysis, R.J.; Funding acquisition, L.H.; Investigation, R.J., L.H. and Z.M.; Methodology, R.J., N.L. and J.-P.G.-E.; Project
administration, R.J. and L.H.; Resources, J.H.; Software, R.J., N.L. and J.-P.G.-E.; Supervision, L.H., Z.M. and J.-P.G.-E.; Validation, R.J.; Visualization, R.J.; Writing—original draft, R.J., L.H. and Z.M.; Writing—review & editing, L.H., Z.M., J.H., N.L. and J.-P.G.-E.

**Funding:** This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic within the National Sustainability Program I (NPU I), grant no. LO1415, and within the Inter-Excellence program (grant no. LTAUSA17175) and partly by the European Space Agency (contract no. 4000115955/15/NL/AF). DART simulations were carried out in part using computational resources supplied by the Ministry of Education, Youth and Sports of the Czech Republic under the Projects CESNET (project no. LM2015042) and CERIT-Scientific Cloud (project no. LM2015085) provided within the program Projects of Large Research, Development and Innovations Infrastructures. The contribution of Zbyněk Malenovský was supported by the Australian Research Council Future Fellowship: *Bridging Scales in Remote Sensing of Vegetation Stress* (FT160100477).

**Acknowledgments:** The authors would like to thank Petr Sloup and Tomáš Rebok for providing code for 3D tree wooden skeleton reconstruction and Jan Novotný for advice regarding terrestrial laser scanning data processing.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Abbreviations

The following abbreviations are used in this manuscript:

<table><tr><td>ASD</td><td>Analytical Spherical Devices</td></tr><tr><td>CASI</td><td>Compact Airborne Spectrographic Imager</td></tr><tr><td>CzeCOS</td><td>Czech Carbon Observation System</td></tr><tr><td>DART</td><td>Discrete Anisotropic Radiative Transfer</td></tr><tr><td>DBH</td><td>Diameter at Breast Height</td></tr><tr><td>FD</td><td>First Derivative</td></tr><tr><td>FLIES</td><td>Forest Light Environmental Simulator</td></tr><tr><td>FRT</td><td>Forest Radiative Transfer</td></tr><tr><td>ICOS</td><td>Integrated Carbon Observation System</td></tr><tr><td>LAD</td><td>Leaf Angle Distribution</td></tr><tr><td>LAI</td><td>Leaf Area Index</td></tr><tr><td>LiDAR</td><td>Light Detection And Ranging</td></tr><tr><td>MSI</td><td>MultiSpectral Imager</td></tr><tr><td>QSM</td><td>Quantitative Structure Model</td></tr><tr><td>RAMI</td><td>Radiation transfer Model Intercomparison</td></tr><tr><td>RMSE</td><td>Root Mean Square Error</td></tr><tr><td>RS</td><td>Remote sensing</td></tr><tr><td>RT</td><td>Radiative Transfer</td></tr><tr><td>RTM</td><td>Radiative Transfer Model</td></tr><tr><td>SAIL</td><td>Scattering by Arbitrary Inclined Leaves</td></tr><tr><td>SASI</td><td>Smaller SWIR Airborne Spectrographic Imager</td></tr><tr><td>TLS</td><td>Terrestrial Laser Scanning</td></tr></table>

## Appendix A. Finding Tree Crown Envelope

The process of finding the crown envelope consists of four steps (illustrated in Figure A1):

1. choosing a voxel $v_c$: the algorithm iterates through every voxel containing at least one shoot position;

2. finding a set of all neighbours of voxel $N_{26}(v_c)$ (26 voxels around the current one);

3. creation of the tree envelope – checking if at least one voxel of the $N_{26}(v_c)$ does not contain a shoot positions; if yes, add the $v_c$ to the set of edge voxels and continue with the next voxel, otherwise continue directly with the next voxel;

4. filling the envelope of a given crown—reduce convexity of the envelope (for more details, see explanation below).
<figure><img src="image_15.png" /><figcaption>**Figure A1.** The process of finding crown envelope. The corresponding steps are described above 1.– 4.</figcaption></figure>
The crown envelope convexity reduction prevents new shoots from being placed inside the crown as they naturally grow at the periphery. This may happen when the air gap size is larger than the voxel size. The following steps were implemented to reduce the inappropriate crown envelope convexity (illustrated in Figure A2):

4.1. finding voxels from the crown envelope in the current horizontal layer—V<sub>c</sub>;

4.2. separating V<sub>c</sub> into lines by x coordinate;

4.3. finding cubes with minimum and maximum y coordinate values—v<sub>min</sub>, v<sub>max</sub>;

4.4. finding gaps between v<sub>min</sub> and v<sub>max</sub>;

4.5. if there are any gaps, fill them with new voxels;

4.6. creating a new set V<sub>cf</sub> and adding new voxels together with voxels from V<sub>c</sub>;

4.7. repeating steps 4.2.–4.5. for the y-axis;

4.8. adding new voxels into V<sub>cf</sub>;

4.9. separating V<sub>cf</sub> in lines by x coordinate;

4.10. finding voxels with the minimum and maximum y coordinate values—v<sub>min</sub>, v<sub>max</sub>;

4.11. adding voxels v<sub>min</sub> and v<sub>max</sub> into V<sub>c new</sub>;

4.12. repeating steps 4.9.–4.11. for the y-axis;

4.13. setting of voxels V<sub>c new</sub> as the new set of envelope voxels for the current horizontal layer.
<figure><img src="image_16.png" /><figcaption>**Figure A2.** The process of crown envelope convexity corrections. This process describes in detail Step 4, consisting in finding the crown envelope. Corresponding steps 4.1. to 4.13 are presented above.</figcaption></figure>
## Appendix B. Assigning Shoots into Age Categories

The tree crown is divided into several horizontal layers ($l_n$) 0.8 m high (approximately one branch whorl). The percentage of the current-year needle category in each layer was assigned according to measurements published by Barták [54] (see Table 1). Differentiation into needle age categories is illustrated in Figure A3 and described by following steps:

1. finding corresponding percentage of current-year shoots in a horizontal layer $l_n$;

2. finding envelope voxels of the tree crown in $l_n$;

3. finding shoot positions inside the envelope voxels;

4. calculating number of shoot positions found;

5. adding shoot positions found into the set of current-year shoot positions;

6. subtracting the shoot positions found from the remaining shoot positions;

7. verifying the percentage of shoot positions in the set of current-year shoot positions; testing if the percentage is higher or lower than the corresponding percentage in $l_n$;

8. deciding to continue or stop the algorithm; if there are not enough current shoot positions, the algorithm continues to the next step, 9; otherwise, it continues to the next layer and repeats steps 1.–7.;

9. extending amount of voxels in the envelope of tree crown towards the tree trunk in $l_n$ is done for x- and y-axes separately (see Figure A4);

10. returning to step 3., then continuing until the end.
<figure><img src="image_17.png" /><figcaption>**Figure A3.** The process of separating and finding crown envelope. The process iterates through every voxel, one by one, so long as there are still points from foliage point cloud.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure A4.** Extension of a crown edge towards the tree trunk. The extension is conducted in two consecutive steps, namely in the x and the y directions. Light green represents voxels contained in tree crown envelope, and dark green represents inside of the crown.</figcaption></figure>

## Appendix C. Angular Distribution of Shoots

The position in which the shoot was supposed to be translated is noted as $L = (l_x, l_y, l_z)$. A shoot was first rotated around the y and z axis. The elevation angles $\beta$, listed in Table 1, are dependent on height ($l_z$). The calculation of the azimuth angle $\gamma$ is based on $l_x$ and $l_y$. It is randomized by an angle $\psi$ getting values within the range $(-\frac{\pi}{4}, \frac{\pi}{4})$ (Equation (A1)).

$$ \gamma = \begin{cases} \arctan(\frac{l_y}{l_z}) + \psi & \text{if } l_x > 0 \\ \arctan(\frac{l_y}{l_z}) + \pi + \psi & \text{if } l_x < 0 \end{cases} \qquad (A1) $$

The vector for shoot translation was calculated as $\mathbf{v} = (L - 0)$, where 0 is the origin of the coordinate system. The transformation matrix was computed as:

$$ \begin{pmatrix} \cos\gamma \cos\beta & -\sin\gamma & -\cos\gamma \sin\beta & v_x \\ \sin\gamma \cos\beta & \cos\gamma & -\sin\gamma \sin\beta & v_y \\ \sin\beta & 0 & \cos\beta & v_z \\ 0 & 0 & 0 & 1 \end{pmatrix} \qquad (A2) $$

Equation (A2) represents the general matrix for shoot models transformation, which was applied to distribute shoots of the two age categories within the reconstructed 3D model of a Norway spruce tree.
## Appendix D. Parametrization of the Basic Turbid DART Tree Representation

Predefined tree shapes available in DART (Fugure 1a) allow constructing spruce trees as composed ellipsoids with precisely defined tree crown dimensions (Table A1). Each tree crown was further divided into three vertical levels. Within each level, different leaf optical properties were defined, tree stem diameter adjusted, and percentage of empty leaf cells (gaps) controlled. Gaps were distributed such that there were more gaps along the stem and at the bottom (which are shaded parts) and thus to mimic natural tree defoliation (Table A2). Leaf optical properties were computed as weighted means of leaf spectra obtained for needles of different categories (sun-exposed needles from the tree top vs. shade-adapted needles from the tree bottom and needles of different age categories ranging from the current year to needles up to 2 years old). Weighting factors (Table A3) were derived from data on distribution of needle age classes within spruce crowns obtained by field measurement conducted at the Bílý Kříž site by Pokorný and Marek [71]. Average leaf optical properties for each crown level are shown in Figure 4a. Leaf angle distribution was in this case fixed and defined as a spherical distribution. Forest scene leaf area index was fully controlled by DART and was set up to be equal to 8.5 m² m⁻².

<figure><img src="image_19.png" /><figcaption>**Table A1.** Basic tree parameters used to construct spruce trees in DART using predefined tree shapes. The figure on the top illustrates the individual tree parameters (adapted from DART manual [30]).</figcaption></figure>

<table><thead><tr><th>Parameter</th><th>Tree 1</th><th>Tree 2</th><th>Tree 3</th><th>Tree 4</th></tr></thead><tbody><tr><td>(1) trunk height below crown (m)</td><td>1.55</td><td>4.6</td><td>0</td><td>2.5</td></tr><tr><td>(2) trunk height within crown (m)</td><td>13.46</td><td>10.41</td><td>16.04</td><td>11.5</td></tr><tr><td>(3) trunk diameter below crown (m)</td><td>0.26</td><td>0.25</td><td>0.27</td><td>0.3</td></tr><tr><td>(4) crown height (m)</td><td>13.45</td><td>10.52</td><td>16.26</td><td>11.65</td></tr><tr><td>(5) intermediate height (m)</td><td>9</td><td>10</td><td>9.5</td><td>9</td></tr><tr><td>(6) crown projection - first axis (m)</td><td>2.58</td><td>2.74</td><td>3.95</td><td>3.5</td></tr><tr><td>(7) crown projection - second axis (m)</td><td>2.38</td><td>2.88</td><td>3.79</td><td>3.3</td></tr></tbody></table>
<table><caption>**Table A2.** Definition of tree crown levels with inner defoliation zone that is illustrated by parameters *a* and *b* (the horizontal axis illustrates tree crown diameter, where 0 represents the tree stem and 1 represents the crown's periphery, adapted from DART manual [30]).</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="3">proportion of full cells</th></tr><tr><th>Top</th><th>Middle</th><th>Bottom</th></tr></thead><tbody><tr><td>level height relative to crown length</td><td>0.17</td><td>0.3</td><td>0.53</td></tr><tr><td>trunk diameter relative to one below</td><td>0</td><td>0.375</td><td>0.75</td></tr><tr><td>percentage of full leaf cells</td><td>0.85</td><td>0.50</td><td>0.30</td></tr><tr><td>start of the inner defoliation zone (a)</td><td>0</td><td>0.3</td><td>0.75</td></tr><tr><td>end of the inner defoliation zone (b)</td><td>1</td><td>1</td><td>1</td></tr></tbody></table>

<table><caption>**Table A3.** Weighting factors used to compute average leaf optical properties for crown vertical levels.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Needle Age</th></tr><tr><th>C</th><th>C+1</th><th>C+2</th></tr></thead><tbody><tr><td>top</td><td>0.70</td><td>0.20</td><td>0.10</td></tr><tr><td>middle</td><td>0.50</td><td>0.25</td><td>0.25</td></tr><tr><td>bottom</td><td>0.10</td><td>0.15</td><td>0.75</td></tr></tbody></table>

## References

1. Mittermeier, R.A.; Myers, N.; Thomsen, J.B.; Fonseca, G.A.B.D.; Olivieri, S. Biodiversity Hotspots and Major Tropical Wilderness Areas: Approaches to Setting Conservation Priorities. *Conserv. Biol.* **1998**, *12*, 516–520. [CrossRef]

2. Myers, N. The biodiversity challenge: Expanded hot-spots analysis. *Environmentalist* **1990**, *10*, 243–256. [CrossRef] [PubMed]

3. Pan, Y.; Birdsey, R.A.; Fang, J.; Houghton, R.; Kauppi, P.E.; Kurz, W.A.; Phillips, O.L.; Shvidenko, A.; Lewis, S.L.; Canadell, J.G.; et al. A Large and Persistent Carbon Sink in the World's Forests. *Science* **2011**, *333*, 988–993. [CrossRef] [PubMed]

4. Schröter, D.; Cramer, W.; Leemans, R.; Prentice, I.C.; Araújo, M.B.; Arnell, N.W.; Bondeau, A.; Bugmann, H.; Carter, T.R.; Gracia, C.A.; et al. Ecosystem Service Supply and Vulnerability to Global Change in Europe. *Science* **2005**, *310*, 1333–1337. [CrossRef] [PubMed]

5. Chazdon, R.L. Beyond Deforestation: Restoring Forests and Ecosystem Services on Degraded Lands. *Science* **2008**, *320*, 1458–1460. [CrossRef] [PubMed]

6. Leuning, R.; Kelliher, F.M.; Pury, D.G.G.D.; Schulze, E.D. Leaf nitrogen, photosynthesis, conductance and transpiration: Scaling from leaves to canopies. *Plant Cell Environ.* **1995**, *18*, 1183–1200. [CrossRef]

7. Trumbore, S.; Brando, P.; Hartmann, H. Forest health and global change. *Science* **2015**, *349*, 814–818. [CrossRef] [PubMed]

8. Chambers, J.Q.; Asner, G.P.; Morton, D.C.; Anderson, L.O.; Saatchi, S.S.; Espírito-Santo, F.D.B.; Palace, M.; Souza, C. Regional ecosystem structure and function: Ecological insights from remote sensing of tropical forests. *Trends Ecol. Evol.* **2007**, *22*, 414–423. [CrossRef] [PubMed]

9. Fassnacht, F.E.; Latifi, H.; Stereńczak, K.; Modzelewska, A.; Lefsky, M.; Waser, L.T.; Straub, C.; Ghosh, A. Review of studies on tree species classification from remotely sensed data. *Remote Sens. Environ.* **2016**, *186*, 64–87. [CrossRef]

10. Hansen, M.C.; Loveland, T.R. A review of large area monitoring of land cover change using Landsat data. *Remote Sens. Environ.* **2012**, *122*, 66–74. [CrossRef]

11. Koch, B. Status and future of laser scanning, synthetic aperture radar and hyperspectral remote sensing data for forest biomass assessment. *ISPRS J. Photogramm. Remote Sens.* **2010**, *65*, 581–590. [CrossRef]
12. Lausch, A.; Erasmi, S.; King, D.; Magdon, P.; Heurich, M.; Lausch, A.; Erasmi, S.; King, D.J.; Magdon, P.; Heurich, M. Understanding Forest Health with Remote Sensing-Part I—A Review of Spectral Traits, Processes and Remote-Sensing Characteristics. *Remote Sens.* **2016**, *8*, 1029. [CrossRef]

13. Lausch, A.; Erasmi, S.; King, D.; Magdon, P.; Heurich, M.; Lausch, A.; Erasmi, S.; King, D.J.; Magdon, P.; Heurich, M. Understanding Forest Health with Remote Sensing-Part II—A Review of Approaches and Data Models. *Remote Sens.* **2017**, *9*, 129. [CrossRef]

14. White, J.C.; Coops, N.C.; Wulder, M.A.; Vastaranta, M.; Hilker, T.; Tompalski, P. Remote Sensing Technologies for Enhancing Forest Inventories: A Review. *Can. J. Remote Sens.* **2016**, *42*, 619–641. [CrossRef]

15. Disney, M.; Lewis, P.; Saich, P. 3D modelling of forest canopy structure for remote sensing simulations in the optical and microwave domains. *Remote Sens. Environ.* **2006**, *100*, 114–132. [CrossRef]

16. Panferov, O.; Knyazikhin, Y.; Myneni, R.B.; Szarzynski, J.; Engwald, S.; Schnitzler, K.G.; Gravenhorst, G. The role of canopy structure in the spectral variation of transmission and absorption of solar radiation in vegetation canopies. *IEEE Trans. Geosci. Remote Sens.* **2001**, *39*, 241–253. [CrossRef]

17. Widlowski, J.L.; Pinty, B.; Gobron, N.; Verstraete, M.M.; Diner, D.J.; Davis, A.B. Canopy Structure Parameters Derived from Multi-Angular Remote Sensing Data for Terrestrial Carbon Studies. *Clim. Chang.* **2004**, *67*, 403–415. [CrossRef]

18. Verhoef, W. Light scattering by leaf layers with application to canopy reflectance modeling: The SAIL model. *Remote Sens. Environ.* **1984**, *16*, 125–141. [CrossRef]

19. Verhoef, W.; Bach, H. Coupled soil-leaf-canopy and atmosphere radiative transfer modeling to simulate hyperspectral multi-angular surface reflectance and TOA radiance data. *Remote Sens. Environ.* **2007**, *109*, 166–182. [CrossRef]

20. Gobron, N.; Pinty, B.; Verstraete, M.M.; Govaerts, Y. A semidiscrete model for the scattering of light by vegetation. *J. Geophys. Res. Atmos.* **1997**, *102*, 9431–9446. [CrossRef]

21. Kallel, A. Vegetation radiative transfer modeling based on virtual flux decomposition. *J. Quant. Spectrosc. Radiat. Transf.* **2010**, *111*, 1389–1405. [CrossRef]

22. North, P. Three-dimensional forest light interaction model using a Monte Carlo method. *IEEE Trans. Geosci. Remote Sens.* **1996**, *34*, 946–956. [CrossRef]

23. Kuusk, A.; Nilson, T. A Directional Multispectral Forest Reflectance Model. *Remote Sens. Environ.* **2000**, *72*, 244–252. [CrossRef]

24. Kobayashi, H.; Iwabuchi, H. A coupled 1-D atmosphere and 3-D canopy radiative transfer model for canopy reflectance, light environment, and photosynthesis simulation in a heterogeneous landscape. *Remote Sens. Environ.* **2008**, *112*, 173–185. [CrossRef]

25. Govaerts, Y.M.; Verstraete, M.M. Raytran: A Monte Carlo ray-tracing model to compute light scattering in three-dimensional heterogeneous media. *IEEE Trans. Geosci. Remote Sens.* **1998**, *36*, 493–505. [CrossRef]

26. Qin, W.; Gerstl, S.A.W. 3-D Scene Modeling of Semidesert Vegetation Cover and its Radiation Regime. *Remote Sens. Environ.* **2000**, *74*, 145–162. [CrossRef]

27. Bailey, B.N.; Overby, M.; Willemse, P.; Pardyjak, E.R.; Mahaffee, W.F.; Stoll, R. A scalable plant-resolving radiative transfer model based on optimized GPU ray tracing. *Agric. For. Meteorol.* **2014**, *198–199*, 192–208. [CrossRef]

28. Gastellu-Etchegorry, J.P.; Demarez, V.; Pinel, V.; Zagolski, F. Modeling radiative transfer in heterogeneous 3-D vegetation canopies. *Remote Sens. Environ.* **1996**, *58*, 131–156. [CrossRef]

29. Gastellu-Etchegorry, J.P.; Lauret, N.; Yin, T.; Landier, L.; Kallel, A.; Malenovský, Z.; Bitar, A.A.; Aval, J.; Benhmida, S.; Qi, J.; et al. DART: Recent Advances in Remote Sensing Data Modeling With Atmosphere, Polarization, and Chlorophyll Fluorescence. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.* **2017**, *10*, 2640–2649. [CrossRef]

30. Gastellu-Etchegorry, J.P. DART User Manual, 2019. Available online: http://www.cesbio.ups-tlse.fr/dart/Public/documentation/contenu/documentation/DART_User_Manual.pdf (accessed on 25 March 2019).

31. Prusinkiewicz, P.; Lindenmayer, A. *The Algorithmic Beauty of Plants*; The Virtual Laboratory, Springer: New York, NY, USA, 1990.

32. Kurth, W.; Sloboda, B. Growth grammars simulating trees–an extension of L-systems incorporating local variables and sensitivity. *Silva Fenn.* **1997**, *31*, 285–295. [CrossRef]
33. Pradal, C.; Fournier, C.; Valduriez, P.; Cohen-Boulakia, S. OpenAlea: Scientific Workflows Combining Data Analysis and Simulation. In Proceedings of the 27th International Conference on Scientific and Statistical Database Management, SSDBM '15, La Jolla, CA, USA, 29 June–1 July 2015; ACM: New York, NY, USA, 2015; pp. 11:1–11:6. [CrossRef]

34. Widlowski, J.L.; Mio, C.; Disney, M.; Adams, J.; Andredakis, I.; Atzberger, C.; Brennan, J.; Busetto, L.; Chelle, M.; Ceccherini, G.; et al. The fourth phase of the radiative transfer model intercomparison (RAMI) exercise: Actual canopy scenarios and conformity testing. *Remote Sens. Environ.* **2015**, *169*, 418–437. [CrossRef]

35. Gruber, F. Morphology of coniferous trees: Possible effects of soil acidification on the morphology of Norway spruce and silver fir. In *Effects of Acid Rain on Forest Processes*; Wiley-Liss: New York, NY, USA, 1994; pp. 265–324.

36. Van Leeuwen, M.; Nieuwenhuis, M. Retrieval of forest structural parameters using LiDAR remote sensing. *Eur. J. For. Res.* **2010**, *129*, 749–770. [CrossRef]

37. Dassot, M.; Constant, T.; Fournier, M. The use of terrestrial LiDAR technology in forest science: application fields, benefits and challenges. *Ann. For. Sci.* **2011**, *68*, 959–974. [CrossRef]

38. Disney, M.I.; Boni Vicari, M.; Burt, A.; Calders, K.; Lewis, S.L.; Raumonen, P.; Wilkes, P. Weighing trees with lasers: Advances, challenges and opportunities. *Interface Focus* **2018**, *8*, 20170048. [CrossRef] [PubMed]

39. Raumonen, P.; Kaasalainen, M.; Åkerblom, M.; Kaasalainen, S.; Kaartinen, H.; Vastaranta, M.; Holopainen, M.; Disney, M.; Lewis, P. Fast Automatic Precision Tree Models from Terrestrial Laser Scanner Data. *Remote Sens.* **2013**, *5*, 491–520. [CrossRef]

40. Hackenberg, J.; Spiecker, H.; Calders, K.; Disney, M.; Raumonen, P. SimpleTree —An Efficient Open Source Tool to Build Tree Models from TLS Clouds. *Forests* **2015**, *6*, 4245–4294. [CrossRef]

41. Sloup, P. Automatic Tree Reconstruction from Its Laser Scan. Master's Thesis, Masaryk University Faculty of Informatics, Brno-Královo Pole-Ponava, Czechia, 2013.

42. Ma, L.; Zheng, G.; Eitel, J.U.H.; Moskal, L.M.; He, W.; Huang, H. Improved Salient Feature-Based Approach for Automatically Separating Photosynthetic and Nonphotosynthetic Components Within Terrestrial Lidar Point Cloud Data of Forest Canopies. *IEEE Trans. Geosci. Remote Sens.* **2016**, *54*, 679–696. [CrossRef]

43. Douglas, E.S.; Martel, J.; Li, Z.; Howe, G.; Hewawasam, K.; Marshall, R.A.; Schaaf, C.L.; Cook, T.A.; Newnham, G.J.; Strahler, A.; et al. Finding Leaves in the Forest: The Dual-Wavelength Echidna Lidar. *IEEE Geosci. Remote Sens. Lett.* **2015**, *12*, 776–780. [CrossRef]

44. Zhu, X.; Skidmore, A.K.; Darvishzadeh, R.; Niemann, K.O.; Liu, J.; Shi, Y.; Wang, T. Foliar and woody materials discriminated using terrestrial LiDAR in a mixed natural forest. *Int. J. Appl. Earth Obs. Geoinf.* **2018**, *64*, 43–50. [CrossRef]

45. Xie, D.; Wang, X.; Qi, J.; Chen, Y.; Mu, X.; Zhang, W.; Yan, G. Reconstruction of Single Tree with Leaves Based on Terrestrial LiDAR Point Cloud Data. *Remote Sens.* **2018**, *10*, 686. [CrossRef]

46. Åkerblom, M.; Raumonen, P.; Casella, E.; Disney, M.I.; Danson, F.M.; Gaulton, R.; Schofield, L.A.; Kaasalainen, M. Non-intersecting leaf insertion algorithm for tree structure models. *Interface Focus* **2018**, *8*, 20170045. [CrossRef] [PubMed]

47. Calders, K.; Origo, N.; Burt, A.; Disney, M.; Nightingale, J.; Raumonen, P.; Åkerblom, M.; Malhi, Y.; Lewis, P. Realistic Forest Stand Reconstruction from Terrestrial LiDAR for Radiative Transfer Modelling. *Remote Sens.* **2018**, *10*, 933. [CrossRef]

48. Côté, J.F.; Widlowski, J.L.; Fournier, R.A.; Verstraete, M.M. The structural and radiative consistency of three-dimensional tree reconstructions from terrestrial lidar. *Remote Sens. Environ.* **2009**, *113*, 1067–1081. [CrossRef]

49. Zheng, G.; Moskal, L.M. Leaf Orientation Retrieval From Terrestrial Laser Scanning (TLS) Data. *IEEE Trans. Geosci. Remote Sens.* **2012**, *50*, 3970–3979. [CrossRef]

50. Bailey, B.N.; Mahaffee, W.F. Rapid measurement of the three-dimensional distribution of leaf orientation and the leaf angle probability density function using terrestrial LiDAR scanning. *Remote Sens. Environ.* **2017**, *194*, 63–76. [CrossRef]

51. Krupková, L.; Marková, I.; Havránková, K.; Pokorný, R.; Urban, O.; Šigut, L.; Pavelka, M.; Cienciala, E.; Marek, M.V. Comparison of different approaches of radiation use efficiency of biomass formation estimation in Mountain Norway spruce. *Trees* **2017**, *31*, 325–337. [CrossRef]
52. Urban, O.; Hrstka, M.; Zitová, M.; Holišová, P.; Šprtová, M.; Klem, K.; Calfapietra, C.; De Angelis, P.; Marek, M.V. Effect of season, needle age and elevated CO2 concentration on photosynthesis and Rubisco acclimation in Picea abies. *Plant Physiol. Biochem.* **2012**, *58*, 135–141. [CrossRef] [PubMed]

53. Homolová, L.; Janoutová, R.; Lukeš, P.; Hanuš, J.; Novotný, J.; Brovkina, O.; Loayza Fernandez, R.R. In situ data supporting remote sensing estimation of spruce forest parameters at the ecosystem station Bílý Kříž. *Beskydy* **2017**, *10*, 75–86. [CrossRef]

54. Barták, M. Struktura Koruny Smrku Ztepilého ve Vztahu k Produkci. Ph.D. Thesis, Kandidátská Dizertační Práce, Ústav Systematické a Ekologické Biologie, Brno, Czechia, 1992.

55. Barták, M.; Dvořák, V.; Hudcová, L. Rozložení biomasy jehlic v korunové vrstvě smrkového porostu. *Lesnictví-Forestry* **1993**, *39*, 273–281.

56. Mesarch, M.A.; Walter-Shea, E.A.; Asner, G.P.; Middleton, E.M.; Chan, S.S. A Revised Measurement Methodology for Conifer Needles Spectral Optical Properties: Evaluating the Influence of Gaps between Elements. *Remote Sens. Environ.* **1999**, *68*, 177–192. [CrossRef]

57. Malenovský, Z.; Albrechtová, J.; Lhotáková, Z.; Zurita-Milla, R.; Clevers, J.G.P.W.; Schaepman, M.E.; Cudlín, P. Applicability of the PROSPECT model for Norway spruce needles. *Int. J. Remote Sens.* **2006**, *27*, 5315–5340. [CrossRef]

58. Yáñez-Rausell, L.; Malenovský, Z.; Clevers, J.G.P.W.; Schaepman, M.E. Minimizing Measurement Uncertainties of Coniferous Needle-Leaf Optical Properties. Part II: Experimental Setup and Error Analysis. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.* **2014**, *7*, 406–420. [CrossRef]

59. Hanuš, J.; Fabiánek, T.; Fajmon, L. Potential of Airborne Imaging Spectroscopy at Czechglobe. *ISPRS 2016*, *XLI-B1*, 15–17. [CrossRef]

60. Richter, R.; Schläpfer, D. Geo-atmospheric processing of airborne imaging spectrometry data. Part 2: Atmospheric/topographic correction. *Int. J. Remote Sens.* **2002**, *23*, 2631–2649. [CrossRef]

61. Verroust, A.; Lazarus, F. Extracting skeletal curves from 3D scattered data. In Proceedings of the Shape Modeling International'99, International Conference on Shape Modeling and Applications, Aizu-Wakamatsu, Japan, 1–4 March 1999; IEEE: Aizu-Wakamatsu, Japan, 1999; pp. 194–201. [CrossRef]

62. Berk, A.; Conforti, P.; Kennett, R.; Perkins, T.; Hawes, F.; van den Bosch, J. MODTRAN6: A major upgrade of the MODTRAN radiative transfer code. In Proceedings of the Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XX, Baltimore, MD, USA, 13 June 2014; International Society for Optics and Photonics: Bellingham, WA, USA, 2014; Volume 9088, p. 90880H. [CrossRef]

63. Malenovský, Z.; Martin, E.; Homolová, L.; Gastellu-Etchegorry, J.P.; Zurita-Milla, R.; Schaepman, M.E.; Pokorný, R.; Clevers, J.G.; Cudlín, P. Influence of woody elements of a Norway spruce canopy on nadir reflectance simulated by the DART model at very high spatial resolution. *Remote Sens. Environ.* **2008**, *112*, 1–18. [CrossRef]

64. Malenovský, Z.; Homolová, L.; Zurita-Milla, R.; Lukeš, P.; Kaplan, V.; Hanuš, J.; Gastellu-Etchegorry, J.P.; Schaepman, M.E. Retrieval of spruce leaf chlorophyll content from airborne image data using continuum removal and radiative transfer. *Remote Sens. Environ.* **2013**, *131*, 85–102. [CrossRef]

65. Yáñez, L.; Homolová, L.; Schaepman, M.E.; Malenovský, Z. Geometrical and structural parametrization of forest canopy radiative transfer by LIDAR measurements. In Proceedings of the 21st ISPRS Congress, Beijing, China, 3–11 July 2008. [CrossRef]

66. Yáñez-Rausell, L.; Malenovský, Z.; Rautiainen, M.; Clevers, J.G.P.W.; Lukeš, P.; Hanuš, J.; Schaepman, M.E. Estimation of Spruce Needle-Leaf Chlorophyll Content Based on DART and PARAS Canopy Reflectance Models. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.* **2015**, *8*, 1534–1544. [CrossRef]

67. Stenberg, P. Simulations of the effects of shoot structure and orientation on vertical gradients in intercepted light by conifer canopies. *Tree Physiol.* **1996**, *16*, 99–108. [CrossRef]

68. Smolander, S.; Stenberg, P. A method to account for shoot scale clumping in coniferous canopy reflectance models. *Remote Sens. Environ.* **2003**, *88*, 363–373. [CrossRef]

69. Rautiainen, M.; Mõttus, M.; Yáñez-Rausell, L.; Homolová, L.; Malenovský, Z.; Schaepman, M.E. A note on upscaling coniferous needle spectra to shoot spectral albedo. *Remote Sens. Environ.* **2012**, *117*, 469–474. [CrossRef]
70. Mõttus, M.; Rautiainen, M.; Schaepman, M.E. Shoot scattering phase function for Scots pine and its effect on canopy reflectance. *Agric. For. Meteorol.* **2012**, *154–155*, 67–74. [CrossRef]

71. Pokorný, R.; Marek, M.V. Test of Accuracy of LAI Estimation by LAI-2000 under Artificially Changed Leaf to Wood Area Proportions. *Biol. Plant.* **2000**, *43*, 537–544. [CrossRef]

<figure><img src="image_20.png" /></figure>

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).