Reconstructing forest canopy from the 3D triangulations of airborne laser scanning point data for the visualization and planning of forested landscapes

ORIGINAL PAPER Reconstructing forest canopy from the 3D triangulations of airborne laser scanning point data for the visualization and planning of forested landscapes Jari Vauhkonen 1,2,3 &Roope Ruotsalainen 1,2 Received: 15 March 2016 /Accepted: 4 November 2016 #The Author(s) 2017. This article is published with open access at Springerlink.com Abstract &Key message We present a data-driven technique to visualize forest landscapes and simulate their future development according to alternative management scenarios. Gentle harvesting intensities were preferred for maintaining scenic values in a test of eliciting public’s preferences based on the simulated landscapes. &Context Visualizations of future forest landscapes according to alternative management scenarios are useful for eliciting stakeholders’preferences on the alternatives. However, conventional computer visualizations require laborious tree-wise measurements or simulators to generate these observations. &Aims We describe and evaluate an alternative approach, in which the visualization is based on reconstructing forest canopy from sparse density, leaf-off airborne laser scanning data. &Methods Computational geometry was employed to generate filtrations, i.e., ordered sets of simplices belonging to the three-dimensional triangulations of the point data. An appropriate degree of filtering was determined by analyzing the topological persistence of the filtrations. The topology was further utilized to simulate changes to canopy biomass, resembling harvests with varying retention levels. Relative priorities of recreational and scenic values of the harvests were estimated based on pairwise comparisons and analytic hierarchy process (AHP). &Results The canopy elements were co-located with the tree stems measured in the field, and the visualizations derived from the entire landscape showed reasonably realistic, despite a low numerical correspondence with plot-level forest attributes. The potential and limitations to improve the proposed parameterization are discussed. &Conclusion Although the criteria to evaluate the landscape visualization and simulation models were not conclusive, the results suggest that forest scenes may be feasibly reconstructed based on data already covering broad areas and readily available for practical applications. Keywords Spatial multicriteria decision analysis .Public participationgeographicinformationsystem(PPGIS) .Remote sensing .Light detection and ranging (LiDAR) .Persistent homology .Alpha shape 1 Introduction Environmental and forestry decision making requires the identification and comparisons of different management Handling Editor: Jean-Michel Leban Contribution of the co-authors Jari Vauhkonen designed the study, implemented all stages related to reconstructing the forest scenes, and drafted the manuscript. Roope Ruotsalainen carried out the survey and analyzed its results with pairwise comparisons and the analytic hierarchy process. He also read and approved the final manuscript. Electronic supplementary material The online version of this article (doi:10.1007/s13595-016-0598-6) contains supplementary material, which is available to authorized users. *Jari Vauhkonen [email protected] Roope Ruotsalainen [email protected] 1 School of Forest Sciences, University of Eastern Finland, Yliopistokatu 7, P.O. Box 111, FI-80101 Joensuu, Finland 2 Department of Forest Sciences, University of Helsinki, Latokartanonkaari 7, P.O. Box 27, FI-00014 Helsinki, Finland 3 Present address: Natural Resources Institute Finland (Luke), Economics and Society, Yliopistokatu 6, P.O. Box 68, FI-80101 Joensuu, Finland Annals of Forest Science (2017) 74:9 DOI 10.1007/s13595-016-0598-6 alternatives based on multiple objectives and stakeholders (Kangas et al. 2008). Integrating information on the impacts of management decisions with preferences of the stakeholders results in a framework called multiple criteria decision analysis (MCDA). When incorporated with a geographic information system (GIS), the resulting applications are nowadays called “spatial MCDA”or “public participation geographic information system”(PPGIS; Sieber 2006), depending on the involvement of the stakeholder. Early forestry applications of spatial MCDA and PPGIS are provided by Store and Kangas (2001) and Kangas and Store (2003), respectively. Projecting future forest landscapes according to alternative management scenarios is a particularly useful MCDA component for participatory forest planning. For example, impacts of growth, management practices, forest (logging) operations, and locations of logging areas may affect the stakeholders’ preferences on the management alternatives. It has long been recognized that even technically simple computer visualizations improve understanding of forest stand dynamics and effects of management decisions (Pukkala and Kellomäki 1988; Burkhart 1992) and facilitate eliciting the related preferences (Tahvanainen et al. 2001; Karjalainen and Tyrväinen 2002). Visualization techniques and applications are reviewed in the forestry context by Mendoza et al. (2006)andFalcão (2008) and techniques for extracting the decision makers’ preferences from the visualizations by Kangas et al. (2008). Recent practical examples of incorporating the information obtained in decision support and GIS frameworks are presented by Warren-Kretzschmar and von Haaren (2014) and Lämås et al. (2015). The visualizations of forested landscapes require balancing between the image realism, the desired resolution, and the measurements available (cf. Bergen et al. 1998;Uusitaloand Orland 2001). Even if extremely photorealistic forest scenes may be obtained by combining botanical models of tree architecture with line graphics or real tree textures presented in virtual reality environments (e.g., Aono and Kunii 1984; Honjo and Lim 2001;Fujisakietal.2008),thetreeand branch-level measurements required by these techniques are time-consuming and expensive. Unless tree-wise inventory data are available, the landscapes need to be populated with trees based on less detailed inventory data. A typical approach is to simulate the locations and appearances of the individual trees based on mean diameter, height, and structure information (e.g., Karjalainen and Tyrväinen 2002). Particularly, the use of mean-based attributes and thus the requirement to simulate the trees will result in generalizations and restrictions of reality (Uusitalo and Orland 2001;Wangetal.2006). Techniques based on remote sensing provide significant advances over conventional forest measurements. Particularly due to its capability to present forest structure as a point cloud mapped in 3D, airborne laser scanning (ALS; also referred to in some instances as “airborne scanning LiDAR”) has become an increasingly popular technique for various forestry applications (Maltamo et al. 2014). However, the use of ALS for visualizing forestry decision making has been rare, to date, even though proposed already by Hill and Veitch (2002), McGaughey and Carson (2003), and Ahlberg et al. (2004). Despite the availability of feature-rich software packages for visualizing ALS point clouds, rendering realistic tree geometry from the point data poses a problem (Simons et al. 2014). Typical approaches to parameterize a forest scene in ALS-based analyses include populating a list of tree locations and dimensions using simple artificial turbid media such as cones, ellipses, voxels (Schneider et al. 2014), or other types of 3D primitives (Koch et al. 2014). A common problem of these approaches is, however, the requirement for a very high data density (e.g., tens of pulses per m 2 ), whereas practically available data sets are typically collected for other purposes such as ground elevation modeling with densities as low as <1 pulses per m 2 (e.g., Nord-Larsen and Schumacher 2012; Villikka et al. 2012). Triangulating point data and subsequent filtering of the triangulations (e.g., Delfinado and Edelsbrunner 1995) has been proposed as an alternative means to represent the 3D canopy surface (Vauhkonen et al. 2014). As opposed to approaches that use fixed image elements and therefore require specifying an artificial pixel or voxel resolution, the triangulations are entirely based on the properties of the point data. Applying filtrations, one can adjust the level of detail in the triangulated point cloud and thus account for canopy gaps and detailed properties existing in the data (see Sect. 2.2.for details). The optimal degree of filtering is found to be forest structure specific (Vauhkonen et al. 2016), whereas less detailed information could suffice for generating realistic visualizations of tessellated landscapes (Vauhkonen 2015). We additionally propose that the topology based on the filtrations could be utilized as a mechanism to simulate changes in the forested landscapes visualized (Sect. 2.3.), which is complicated based on simulated forest stand dynamics (cf. Dreyfus 2012). Even though the low resolution obtainable by the practically available sparse data most likely restricts the detail of the visualizations, the decision maker may be expected to benefit from representing real tree geometry (cf. Uusitalo and Orland 2001) and visual descriptions (cf. Tahvanainen et al. 2001) instead of simulated and verbally or numerically described forest scenes. Another practical benefit of our approach is that required data can be feasibly extracted from any location covered by an ALS campaign such as those designed for collecting sparse density data for ground elevation modeling. Although the suitability of such data for forest inventories was verified (Nord-Larsen and Schumacher 2012; Villikka et al. 2012), the benefits of using the spatially explicit 3D scenes obtainable have not been fully employed despite the potential identified earlier (Hill and Veitch 2002; McGaughey and 9 Page 2 of 13 Annals of Forest Science (2017) 74:9 Carson 2003; Ahlberg et al. 2004). Therefore, rather than developing another virtual reality platform aiming at improved photorealism with high costs of collecting the required data, we present a wireframe model on how to extract information from the visualized forest scenes for a spatial MCDA to open up discussion on the potential and limitations of the broadly available sparse ALS point data to support forestry decision making. The purpose of this study is to describe a technique for visualizing forest scenes from sparse density ALS data and evaluate their potential for forestry decision support, making a difference between (i) close-range and (ii) landscape levels. The close-range visualizations correspond to a within-forest view and up to a plot scale, whereas the landscape level shows forest areas in scales of tens to hundreds of hectares corresponding to typical properties managed by a single forest owner. The forest scenes were obtained applying a triangulating and filtering technique to ALS data acquired originally for ground elevation modeling. Regarding (i), we assessed the correspondence of the forest canopy elements modeled as tetrahedra and obtained by filtering the triangulations with forest attributes observed in the field. Regarding (ii), we introduce an intuitive concept based on using the filtrations to simulate changes to canopy biomass. The initial and simulated landscapes are visualized and demonstrated in eliciting public’s preferences on scenic values resulting from biomass harvests applied with varying intensities. 2 Material and methods 2.1 Study area and data The study area is a typical, managed boreal forest in eastern Finland. An area of 1000 × 1200 m located approximately on 62° 31.5′N, 30° 11′E was selected for the visualizations due to the presence of a forest and lake mosaic, which likely produced high recreational and scenic values over the area. The data for the visualizations included the locations of the lakes and ALS data extracted from a topographic database. In addition, the ALS data were extracted for a number of field sample plots existing in the area, and the field measurements colocated with the ALS data were used to assess the correspondence of the modeled canopy with respect to the field data. The plots were located up to 2.5 km apart the visualized area and were earlier used to study the spatial pattern and diameter distribution of the trees from both the ecological and economic perspectives. Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies [L.] Karst.), with a minor proportion of deciduous species, were present in the plots, whereas the area used for the visualizations was assumed to be nearly purely pine-dominated. The ALS data used in the study were acquired by the National Land Survey of Finland as a part of their data acquisition campaign for creating a nationwide ground elevation model for Finland. The data were downloaded from a file service (https://tiedostopalvelu.maanmittauslaitos. fi/tp/kartta?lang=en), from which they are available for free and with extensive permissions of use. The data were acquired on April 30, 2012, with Leica ALS60 scanner operated in a multipulse mode. The flying altitude was 2350 m, yielding a nominal pulse density of 0.8 m −2 . The data provider had detected and classified the ground level, on which the normalization of the vegetation height values was based. As the data are meant specifically for ground elevation modeling, we assumed the accuracy of this classification to be appropriate for our purposes. The analyses were focused only on the first echoes (i.e., “only”and “first of many”of up to four echoes recorded per pulse), aiming to obtain the main information from the data (Vauhkonen et al. 2014), while retaining most generalization abilities over sensors that record a different number of echo categories (e.g., Næsset 2014). The field plots were measured in May–June, 2010. All trees with either diameter at breast height (DBH) ≥4cmor height ≥4 m were mapped for locations and measured for species, DBH, and height. The measurements are described in more detail by Packalén et al. (2013). For simplicity, the plot size used in this study was standardized to 400 m 2 ,rather than varying from 400 to 900 m 2 as in the earlier publications. The standardization was achieved by creating rectangular windows of 20 × 20 m with the center and orientation corresponding to the original plot. The trees located within this window based on the measured tree coordinates were extracted for the analyses. Plot-level basal area was computed based on summing from the diameter measurements. Dominant height was determined as the mean height of 100 thickest trees per hectare (four trees per plot). Individual stem volumes were estimated by models of Laasasenaho (1982), employing the DBH, height, and tree species as predictors. The models for birch were used for all deciduous trees. The spatial pattern of the trees was assessed by means of the Clark-Evans index (CEI; Clark and Evans 1954) of the aggregation of a point pattern. The CEI values were computed based on the tree coordinates using the spatstat package of the R statistical computing environment (Baddeley and Turner 2005), applying an edge correction proposed by Donnelly (1978). General descriptive statistics of the field data are shown in Table 1. 2.2 Modeling the initial forest canopy using triangulations In our approach, “visualization”essentially refers to visualizing forest canopy geometry based on the 3D measurements collected by an ALS sensor (e.g., Sun and Ranson 2000). The geometry is composed of n-dimensional basic elements (points, Annals of Forest Science (2017) 74:9 Page 3 of 13 9 edges, facets, tetrahedra), which are referred to using generic term “simplex”unless an element of a specific dimension is particularly denoted. The main canopy elements are tetrahedra obtained by subdividing the underlying space of the 3D measurements, which is called triangulating. The entire triangulation, referred to below also as “complex,”thus represents both the canopy and empty space, which need to be divided to separate subcomplexes constructing a filter. The triangulations were filtered using a well-known computational geometric or topological concept called 3D alpha shapes (Edelsbrunner and Mücke 1994), while an appropriate degree of filtering was determined by analyzing persistence of the observed topology, i.e., topological persistence or persistent homology (Edelsbrunner et al. 2002). All computations were based on Delaunay triangulations and implemented with C++ and open-source libraries (Morozov 2012;Daetal.2013;Pion and Teillaud 2013). The applied workflow is described below, but the mathematical formalism is intentionally left to the original publications. Illustrations of all the concepts applied to ALS data are presented by Vauhkonen (2015), for example. Following earlier analyses with ALS point data (Vauhkonen et al. 2014,2016), the landscape was tessellated to a regular grid with a cell size corresponding to the plot size (20 × 20 m). The ALS data of each cell or plot was triangulated and filtered separately. The problem related to using the alpha shapes is the determination of a proper value for alpha (α), which is a threshold for the squared radius of the circumscribing sphere of each simplex and therefore a criterion determining which simplices of the full triangulation belong to its current subcomplex, i.e., α-complex or α-shape. Here, the value of αwas selected based on filtration obtained by ordering the α-complexes according to the values of α(see also Delfinado and Edelsbrunner 1995). During the filtration, i.e., by increasing the value of α, simplices of a dimension djoin together with other d-dimensional simplices to form structures with a higher dimensionality. Particularly, when the value of αallows two points (d=0) to connect, an edge of d= 1 between these two points is formed. The edges further join to form facets (d= 2) and tetrahedra (d= 3), the latter being the highest dimension considered here. Each simplex can thus be described by its persistence in the structure formed. This index of persistence was obtained as the absolute difference between the index values of αcausing the birth and death of simplices of the given dimension. Following the logic of the previous paragraph, a birth and death diagram was computed for indices with d= 1 and d=2. Only these dimensions need to be considered, since those with d=0ord= 3 do not born or die, respectively, in this process. In the diagram, the most interesting observations are those offdiagonal: These simplices are interpret to cause the most fundamental changes to the triangulated structure, typically at the lowest values of α, while the diagonality of the observations indicates stability (persistence) of those features. The maximum off-diagonal αat the death of both d= 1 and d= 2 was hypothesized to reflect the primary persisting structures of the obtained triangulation and was thus used as the αvalue determining the degree of filtering and thus the initial state of the forest canopy. The rationality of the obtained canopy was evaluated using the field data. The observed tree (stem) positions and dimensions were visualized with the canopy to assess the co-location of these elements. Descriptive characteristics were extracted from the triangulation-based canopy models and related with the forest attributes measured in the field (Table 1). The value of the parameter αresulting from the persistent homology analysis and the total tetrahedral volume, percent of the plot area covered by the tetrahedral, and top height of the underlying space of the specified α-complex were extracted, assuming these to be related to the total tree attributes of a field plot. Additionally, the number of connected tetrahedra with at least one joint edge was computed, assuming this to be related to the number of tree patches within the plot. The correspondence was assessed by means of the coefficient of determination (R 2 ). 2.3 Simulating changes to the canopy and assessing their impacts to the scenery 2.3.1 Visualizing the landscape and harvests of canopy biomass The full landscapes were visualized using triangulation representation (TriRep) and triangular surface plot (trisurf)functions implemented in Matlab, version R2012a (MathWorks, Table 1 Central characteristics of the 71 field plots Attribute Mean Std Min Max Stem volume, m 3 /ha 199.6 73.2 69.6 482.2 Basal area, m 2 /ha 24.9 6.6 12.2 43.7 Number of stems, 1/ha 1275.3 590.8 400.0 2875.0 Basal-area weighted mean diameter, cm 21.6 5.2 12.8 32.3 Dominant height, m 20.4 4.1 13.5 32.5 CEI, index value 1.09 0.16 0.69 1.41 Std standard deviation, Min minimum, Max maximum 9 Page 4 of 13 Annals of Forest Science (2017) 74:9 Inc., Natick, MA, USA). The landscapes were composed of the adjacent 20 × 20 m cells processed according to Sect. 2.2. Due to depicting solely canopy geometry, certain simplifications to the visualizing other elements were adopted. The number or the spatial pattern of the tree stems could not be precisely modeled, which is a typical result of ALS-based analyses (e.g., Packalén et al. 2013). However, the landscape-level visualizations were found to be more realistic with at least some shadowing caused by the stems. A number of the highest tetrahedra per cell were thus accompanied with simulated tree stems. The number of stems was determined as ceil(n p ×f×V p /H p ), where n p was the number of tetrahedra, V p was the total volume of these tetrahedra, H p was the maximum ALS height value of the cell p,andfwas a unit area factor, in this case 1/400. The diameters of the stems (in meters) were simulated as H t /50, where H t is the maximum height of the considered stem tin meters. The applied rules resulted in too few and too large trees, in general, but provided the landscapes with such a realism that only the largest tree stems were visible when looking from a distance. The produced landscapes were colored as pure pine forests, i.e., species recognition was not attempted due to the low proportion of other species. The ground level was depicted as perfectly flat, since the fluctuations in the terrain elevation were known to be minor in this area. In addition to visualizing the initial canopy, reductions to the canopy volume were simulated to portray harvests with varying intensities. In the approach, simplices belonging to the canopy were re-assigned as empty space by adjusting the value of αparameter (see Sect. 2.2.) downwards. As a result, the initial α-complexes were subtracted starting from the largest tetrahedra, assuming this to correspond with harvests with varying retention levels based on removing the largest trees from the forest. The harvests were focused on one compartment with a dense initial canopy. Altogether, five harvesting intensities were defined by multiplying the initial αby a factor of 0.95, 0.75, 0.68, 0.6, and 0.5, which corresponded to an average reduction of 17, 71, 83, 93, and 99%, respectively, of the initial canopy volume (11,673 m 3 /ha on the area of the full compartment, on average). The visualizations of the full landscape obtained as a result of the harvesting were labeled as images A–E according to the increasing harvesting intensity and are presented as Electronic Supplementary Material (Online Resource 1). 2.3.2 Pairwise comparisons of the simulated landscapes To obtain preference information on how the harvesting affected the landscape, a pairwise comparison of the visualized landscapes was implemented using the SurveyMonkey® online platform (http://www.surveymonkey.com/). In the survey, two simulated landscapes were shown at the time, asking a question of which of the two images was better with respect to the recreational use and the scenic values of the landscape and how much better it was. The respondent was asked to express his/her opinion on whether either image was (i) equally good, (ii) slightly better, (iii) clearly better, (iv) considerably/strongly better, or (v) absolutely better compared to the other. In the following text, this decision is called the priority ratio of the respondent between the two landscapes. A link to the survey was submitted to the mailing list of the forestry students of the University of Eastern Finland and accompanied with brief instructions (translated from Finnish): “The central compartment of [this] landscape is marked to be harvested. A walking trail and a ‘laavu’(a temporary shelter for field lunch or camping) are situated in the area of the compartment. The harvesting intensity is not yet determined, but the aim is to maintain the recreational value of the forest as high as possible. In the survey, you are asked to compare two alternative landscapes and select which harvesting intensity is better suited with respect to the recreational use and scenic values. The compartment to be harvested is delineated with a red boundary in the image below. The other parts of the landscape will not be affected. Start the survey by clicking ‘Next image pair’button.” The survey was held open between 21 January and 12 February, 2015. The survey yielded altogether 57 responses, of which 7 were incomplete (no response to one or more comparisons). Altogether, 50 responses were thus analyzed in the further study. 2.3.3 Preference elicitation On the basis of the pairwise comparisons (Sect. 2.3.2.), relative priorities of the scenic values with respect to the harvest intensities were estimated using the eigenvalue technique of analytic hierarchy process (AHP; Saaty 1977,1980). This technique was used to convert verbal evaluations to ratioscale numerical values in a forestry decision support application bearing a close resemblance to our study (Kangas et al. 1993). In the approach, the priority ratios i–v(Sect. 2.3.2.) are translated into numerical values of 1:1, 3:1, 5:1, 7:1, and 9:1, respectively. The values are arranged into a reciprocal matrix A, and using the matrix as input, the numeric weights or priorities of the pairwise comparisons are computed as the right eigenvector of the largest eigenvalue (λ max )ofA(Saaty 1977). To assess the coherence of the pairwise comparisons, a consistency index (CI; Saaty 1980), which estimates the level of consistency with respect to the entire comparison process, is derived by relating the value of λ max with the number of comparisons made. The observed CI is further compared with Annals of Forest Science (2017) 74:9 Page 5 of 13 9 an average CI obtained by randomly generated comparisons in a matrix with a size corresponding to A(Saaty 1980). The result is a consistency ratio (CR), which measures the coherence on a scale of 0.1. According to a general rule of thumb, comparisons with CR > 0.1 are deemed inconsistent. 3Results Except for top height, the triangulated and filtered canopy models were in low correspondence with any forest attribute when evaluated with respect to the field measurements. Nevertheless, the R 2 values for the total volume and top height of the canopy model in particular were statistically significant (Table 2), and the forest structure affected the models obtained in the way it was expected. The spatial pattern of the trees, quantified in terms of CEI, significantly affected the values of the αparameter (Table 2), indicating that a different αvalue was required for different spatial patterns. The total tetrahedral volume of the reconstructed canopy was significantly related to all other size and structure attributes expect for the stem number (Table 2), and the volume increased according to increasing values of these attributes. The number of connected tetrahedra and the coverage of the tetrahedra relative to the plot area however had a low degree of determination with any of the field attributes. The top height of the canopy model had the best overall correspondence with the field attributes and a strong linear relationship especially with the dominant tree height. Noting that no allometric knowledge was used in the construction of the canopy, the models based on the geometry and topology of the point data still reflect the variation in the forest properties and aggregated stem attributes. Based on a visual assessment, the canopy elements of the coniferous dominant canopy in particular were fairly accurately co-located with the tree stems measured in the field (Fig. 1). The positions of the canopy elements were also somewhat in line with the spatial patterns of the field measured trees, particularly those showing clustered patterns. Nevertheless, Fig. 1also shows the limitations of the proposed technique toward close-range visualizations. The level of detail and the representation of the deciduous and understorey trees are in particular limited, the reasons of which are further discussed in Sect. 4. Despite the low detail in close-range visualizations (Fig. 1), the visualizations derived for the entire landscape showed a reasonable level of detail and realism (Fig. 2). For instance, roads, electric lines, or otherwise open areas could be pointed out from the image. Further, despite the potential to edge effects, i.e., incorrectly cutting the ALS data due to an edge of a cell (Fig. 1), the borders of the 20 × 20 m cells did not show in the visualization of the full landscape (Fig. 2). The canopy biomass reductions simulated to the landscape were found to reasonably resemble effects of real-world harvesting operations (Fig. 2). However, the simulated harvests clearly had a more considerable effect on the horizontal coverage than on the height of the modeled canopy (Fig. 3). Especially, the top height changed very moderately based on harvesting intensities A, B, and C according to Fig. 3, which shows the effects of the simulated reduction of parameter αto the vertical and horizontal structure of the landscape. According to the survey, the landscape treated with heavy harvesting intensities (images E and D) was preferred by none of the responders, while lower intensities (images C, B, and A) were selected as most preferable by 10 (20%), 16 (32%), and 24 (48%) responders, respectively. Due to these preferences, the (accumulation of) relative priority between the landscapes was significantly indifferent for different types of responders (Fig. 4) and the general trend was difficult to model without additional information. According to the AHP, the consistencies of the responses varied considerably. Only 10 (20%) of the pairwise comparisons had a CR < 0.1, while 30 (60%) and 10 (20%) resulted in CRs between 0.1–0.3 and >0.3, respectively. The maximum and mean CR values were 0.99 and 0.22. When examined against the relative priorities (Fig. 4), the inconsistencies were neither clustered in any preference group nor reflected in any other particular way, however. 4Discussion According to the results presented above, the canopy models and visualizations based on the sparse density ALS data are Table 2 The degree of determination (R 2 )between selected field attributes and alpha value (α), total volume (V), number of connected tetrahedra (N), coverage (C), and top height (H) of the canopy models Field observation αVNCH Stem volume, m 3 /ha 0.01 0.26* 0.07* 0.09* 0.54* Basal area, m 2 /ha 0.00 0.11* 0.06* 0.07* 0.16* Number of stems, 1/ha 0.02 0.02 0.00 0.03 0.19 Basal-area weighted mean diameter, cm 0.05 0.27* 0.04 0.05 0.78* Dominant height, m 0.02 0.32* 0.06* 0.07* 0.94* CEI 0.06* 0.01 0.05 0.06* 0.00 Asterisk marks significance of the correlation at the 95% confidence level, while R 2 > 0.5 are printed in italics 9 Page 6 of 13 Annals of Forest Science (2017) 74:9 mainly suited for landscape-level analyses. The low degree of determination between the close-range models and forest attributes indicates that either a field sample to optimize the canopy model (Vauhkonen et al. 2014,2016) or a completely different approach based on higher density data obtained by ALS or terrestrial laser scanning is required for detailed modeling. On the other hand, the proposed approach appears interesting due to the possibility to generate the visualizations from very sparse ALS data that often readily exist and the inherent property of the modeling technique to simulate changes to the landscape by adjusting the filtration parameter α. The obtained visualizations are affected by both input data and proper selection of the value for α, as discussed below. As opposed to conventional forest visualization techniques that require laborious tree-by-tree measurements or the use of simulators to generate these observations, the present study requires a wall-to-wall coverage of ALS data over the landscape to be visualized. While the acquisitions of ALS data for ground elevation modeling nowadays allow applications of these data for vast areas, properties such as the scarcity of the input point data may restrict some forestry analyses (e.g., Vauhkonen et al. 2016). As seen from Fig. 1of this study, for instance, the true canopy height was underestimated due to the penetration of the pulses through the top canopy (Gaveau and Hill 2003). Nevertheless, the top height extracted from the triangulation-based canopy models had a strong relationship with the dominant height (Table 2), which was comparable to the height metrics extracted from the initial point data (detailed results not shown). No information related to the spatial pattern of the trees or tree number could be derived, which is in line with earlier studies based on ALS data (e.g., Packalén et al. 2013). The percent coverage of the tetrahedra could be compared to the canopy cover metric computed from the point data as the proportion of the echoes above a certain height threshold to all echoes (Korhonen et al. 2011). This metric had a correlation of 0.8 with the tetrahedral coverage, but similarly low degrees of determination with the fieldmeasured forest attributes. Since the data were captured under a leaf-off period of the vegetation phenology, less first-return data were obtained from the deciduous than coniferous trees, even if those were located in the dominant canopy (Fig. 1; see also Villikka et al. 2012). Also, only limited information was obtained from the understorey. Although some researchers have detected more trees under leaf-off than leaf-on conditions, when focusing on individual tree detection in dense ALS data (Duncanson et al. Fig. 1 Examples of visualized plots with clustered (CEI = 0.7, above), random (CEI = 1.0, middle), and regular (CEI = 1.4, below) patterns of tree stems. The locations and dimensions of the cylinders correspond with the tree stems measured in the field. The brown and gray cylinders represent scots pine and deciduous trees, respectively. The interval between tick marks corresponds to 5 m in the vertical axis and 2 m in the plane Annals of Forest Science (2017) 74:9 Page 7 of 13 9 2014), also the results of Korpela et al. (2012) indicate that direct measures of the understorey are difficult to obtain due to transmission losses occurring in the upper canopy (see also Maltamo et al. 2005). For the reasons stated above, the spatial patterns of the deciduous and understorey trees may be incorrect in some parts of the landscape, but this was not assumed to be a major problem in the pine-dominated test area considered in our study. Whether attempted in more complex canopies in further studies, however, considerations on, e.g., the species and pulse penetration properties are required. Whether the limitations described in the previous paragraph are accepted, the presented approach is proposed feasible for visualizing landscapes such as that presented in Fig. 2. The methodology presented here allows a description that is based on real canopy elements, as extracted from the triangulations based on ALS data, but not explicit prescription of leaf or shoot locations and orientations, for example. Therefore, regarding the degree of detail in describing the distribution of canopy elements, the proposed approach is a compromise between artificial tree/crown-level turbid media parameterized by a list of tree locations and dimensions (e.g., Schneider et al. 2014) and more realistic 3D models of vegetation elements potentially obtainable by terrestrial laser scanning. An application to forest management planning is described in the present study, but the modeling approach could also be suited for other applications requiring geometrically explicit parameterizations of forest landscapes. An example is light interaction and radiative transfer modeling (e.g., Schneider et al. 2014), in which incorrect assumptions on crown shapes and inappropriate voxel grid resolutions are reported to result in significant errors in retrieved crown parameters (Calders et al. 2013; Widlowski et al. 2014). Visualization is a straightforward application for the triangulations and filtrations of ALS data collected over a forested area, which are earlier tested in various forest inventory applications (Vauhkonen et al. 2014,2016; Vauhkonen 2015). As opposed to the typical “surface modeling”approaches using 2D to 2.5D raster images, the use of triangulations and filtrations provides means of constructing genuine 3D parameterizations of forested landscapes instead of just representing counts of height values within pixel or voxel cells. Besides triangulations, geometric modeling techniques that could be considered for solving the presented problem are an iterative surface wrapping technique based on delineated tree crowns (Kato et al. 2009), Fig. 2 The full landscape visualized with the approach proposed. The sub-figures on the right give examples on harvests simulated to the compartment in the center of the landscape. All five simulated landscapes are presented as Electronic Supplementary Material (Online Resource 1) Fig. 3 The effects of the simulated reduction of parameter αto the top height and proportion of the plot area covered by the tetrahedra of the modeled canopy. The dots and lines indicate the mean and standard deviation, respectively, within the area. Harvesting intensity 0 corresponds to the initial canopy 9 Page 8 of 13 Annals of Forest Science (2017) 74:9 3D clustering followed by convex polytope reconstruction (Gupta et al. 2010), stacking horizontal slices of the height values (Tang et al. 2013), or voxel-based techniques (Schneider et al. 2014; see also Koch et al. 2014). All these methods involve certain parameters to be set by the operator, and the obtained results depend on these. Further, the methods are so far tested with data densities varying from 4 to 5 m −2 (Gupta et al. 2010)to10–20 m −2 (Kato et al. 2009;Tangetal. 2013; Schneider et al. 2014) as opposed to 0.8 m −2 of this study. The data studied here would have resulted in an extremely coarse voxel grid resolution. On the other hand, requiring data with a higher density would mean distinct data acquisitions and therefore considerably higher inventory costs. Applying filtrations, one can adjust the level of detail in the triangulated point cloud and thus account for canopy gaps and detailed properties existing in the data, which is controlled by filtration parameter α. The approaches to determine the filtration parameter have included using the most feasible population-specific fixed value like in Vauhkonen et al. (2012), who evaluated a range of α’s, but rather than selecting one particular value, they used metrics quantifying the difference of the shape obtained with the range of α’s to the convex hull of the point data, which corresponds to α-shapes with α→∞. Vauhkonen et al. (2014,2016) optimized the degree of filtration with respect to field-measured forest attributes and predicted this degree by means of linear regression. In the present study, the selection of αwas based on analyzing the persistent homology of the filtrations (Edelsbrunner et al. 2002). This approach resulted in less accurate parameter estimates (e.g., R 2 = 0.26 of the volumes in the present study vs. 0.49–0.83 reported in the previous studies) but eliminated the need for a field training data required for the optimization (Vauhkonen et al. 2014,2016). The persistent homology approach was tested to replace the calibration field data by analyzing the inherent topological properties of the point data. A simple approach based on the birth and death diagrams of the various simplices (Edelsbrunner et al. 2002) turned out to produce reasonably good information. The persistent homology of vegetation point clouds should be more comprehensively analyzed after this initial test. For example, one could derive the Euler characteristic and the Betti numbers from the filtrations (Robins 2002) and use them in a similar manner than the birth and death diagram was used. One could further attempt to quantify the obtained filtration by other geometric features than total volumes such as edge lengths or facet areas. The approach has a profound mathematical formalism on topological connectivity and implementations based on efficient data structures are available (Pion and Teillaud 2013;Daetal.2013). A definite difference and a strength of the triangulationbased approach over alternative 3D modeling or description techniques is the ability to simulate changes to the geometric models obtained. This is an inherent property of the technique adopted and may be simply implemented by adjusting the degree of filtering, but, to the best of our knowledge, this option has not been utilized elsewhere beyond the present study. Here, it was used to depict impacts of biomass harvesting to the scenic values of the area by simulating harvests with varying intensities. According to the high CR values, some incoherence existed in the pairwise comparisons, which can possibly be explained by a comparison of Figs. 3and 4. According to Fig. 3, the reduction of the αvalue changed the top height of the canopy models based on harvesting intensities A, B, and C moderately, whereas the reduction of the top height was clearly more pronounced for harvests D and E. This corresponds well with the dispersion between A, B, and C as the most preferred harvesting intensity and the clearly lower preference toward D and E (Fig. 4). Although the canopy coverage changed more than height (Fig. 3), the moderate change in the top height may have dominated especially in the Fig. 4 The distribution (left) and accumulation of the relative preferences between harvesting intensities. The thick line depicts the mean of all preferences. The dashing of the individual preferences varies according to which harvesting intensity was preferred most Annals of Forest Science (2017) 74:9 Page 9 of 13 9