VOXEL BASED REPRESENTATION OF FULL-WAVEFORM AIRBORNE LASER SCANNER DATA FOR FORESTRY APPLICATIONS

The advantages of using airborne full-waveform laser scanner data in forest applications, e.g. for the description of the vertical vegetation structure or accurate biomass estimation, have been emphasized in many publications. To exploit the full potential offered by airborne full-waveform laser scanning data, the development of voxel based methods for data analysis is essential. In contrast to existing approaches based on the extraction of discrete 3D points by a Gaussian decomposition, it is very promising to derive the voxel attributes from the digitised waveform directly. For this purpose, the waveform data have to be transferred into a 3D voxel representation. This requires a series of radiometric and geometric transformations of the raw full-waveform laser scanner data. Thus, the paper deals with the geometric aspects and describes a processing chain from the raw waveform data to an attenuationcorrected volumetric forest stand reconstruction. The integration of attenuation-corrected waveform data into the voxel space is realised with an efficient parametric voxel traversal method operating on an octree data structure. The voxel attributes are derived from the amplitudes of the attenuation-corrected waveforms. Additionally, a new 3D filtering approach is presented to eliminate non-object voxel. Applying these methods to real full-waveform laser scanning data, a voxel based representation of a spruce was generated combining three flight strips from different viewing directions. * Corresponding author


INTRODUCTION
Since commercial small-footprint full-waveform laser scanner became available, the suitability of their data is analysed for different applications.Especially in forestry, the data provide a large gain of information.This allows for a more accurate differentiation of the vegetation structure (e.g.inner canopy and low vegetation).Various publications deal with the use of fullwaveform data for forestry applications such as single tree detection, biomass estimation, tree species classification and vegetation structure characterisation.In the early stages, the data are analysed based on discrete 3D points and additional parameters extracted from the waveforms by peak detection methods like the Gaussian decomposition (e.g.Gupta et al., 2010, Heinzel & Koch, 2011).Moreover, approaches exist based on 3D raster representations derived from discrete points (Reitberger et al., 2009, Leiterer et al., 2012, Bienert et al., 2014, Schilling et al., 2012).So far, only a few publications can be found dealing with the analysis of the entire waveform information although it is very promising for the volumetric modelling of the vertical vegetation structure (Lindberg et al., 2012).To exploit the full potential of full-waveform laser scanning data, methods have to be developed to convert the ALS waveform data in a 3D raster representation which is the subject of this paper.Concerning this, one possible approach is the location of the 3D coordinates of the waveform samples in a defined 3D grid (e.g.Wu et al., 2012, Buddenbaum et al., 2013, Hermosilla et al., 2014).In contrast, this paper presents an octree based ray tracing method which allows efficient waveform integration and additionally the provision of information about the laser pulse history.Assuming an infinitely thin ray by neglecting the laser beam divergence, the integration of the waveforms is realised applying the parametric voxel traversal method published by Revelles et al. (2000).To achieve a volumetric representation with high geometric quality, an accurate determination of the ray from the scanner origin along the direction vector of the laser pulse is required.Hence, the paper also specifies the access to the raw full-waveform laser scanning data and the required geometric transformations of the data.Beyond the geometric aspects, radiometric transformations have to be considered.Due to attenuation effects caused by interactions of the laser pulse with the canopy, a correction of the backscattered signal is required to emphasize sample amplitudes of lower reflections.To achieve a higher radiometric quality of the voxel space, the raw waveforms are corrected using the integral method proposed by Richter et al. (2014).Thus, the voxel attributes are derived from the amplitudes of the attenuation-corrected waveforms.We currently use the received waveform instead of the differential backscatter cross section, thus implying the assumption of an infinitely short emitted laser pulse.To consider redundant voxel entries occurring in the overlapping areas of different flight strips, two scan angle-based methods are compared to the maximum amplitude approach proposed by Hermosilla et al. (2014).Noise present in the signal is reduced with a 3D filtering method which operates directly in the voxel space.The results are compared to common approaches (e.g.Persson et al. 2005, Wu et al., 2012, Hermosilla et al., 2014).The methods were implemented and tested on a real full-waveform dataset.The results of the geo-referencing, voxel space generation and filtering are presented and analysed.

DATA
The full-waveform data set comprises four flight strips with about 75 % overlap which were acquired in March 2015 covering parts of the Tharandt Forest.Additionally, a fifth strip was flown which is located similar to the second strip.The data were recorded with a Riegl LMS Q680i from 600 m height above ground with 400 kHz pulse repetition rate, ±30° field of view and a laser beam divergence of 0.5 mrad.This yields to a laser footprint size of about 30 cm in nadir.The waveforms are discretised with a sampling interval of 1 ns corresponding to approximately 15 cm sampling distance.The raw full-waveform laser scanning data of each flight strip are separately stored in a binary SDF file format.The respective trajectory data providing time-solved GPS positions and measurements of the inertial measurement unit (IMU) are offered in a binary POF file.

DATA PREPARATION
This section describes the steps from single waveforms consisting of time stamps and amplitudes to geo-referenced waveforms which can be located globally and subsequently transferred into a voxel space representation.All algorithms are implemented in C++ using the RiWaveLIB provided by RIEGL for accessing the laser scanner data.

Data reading
The raw full-waveform laser scanning data are stored in binary SDF files.Because there is less information about the file format, Jalobeanu & Gonçalves, (2012) present an approach for reconstructing the file content.Figure 1 shows the generalised structure of the SDF file format.The file header contains general information like the applied laser scanner type and its serial number, the group velocity of an emitted pulse (v) and the sampling interval of the waveforms (t interval ).The scan lines are separated by control records which contain housekeeping data and available information about GPS synchronisation.Each emitted laser pulse can be retrieved as pulse record in the SDF file.A pulse record (Fig. 2) comprises header information including origin (o SOCS ) and direction vector (d SOCS ) of pulse emission, the number of recorded sample blocks, and the instant of time when the data sampling was initiated as internal time (time of start of range gate t sorg ) and GPS time (t sorg GPS ).The discretised backscattered waveforms are stored in sampling blocks.Each sampling block contains e.g.information about the number of samples (n) and the time of the first recorded sample of the block (starting time of sampling block t sosbl ).In general, the first sampling block represents the emitted waveform.The reflected waveform can be represented by multiple sampling blocks if the height of the reflecting object exceeds the fixed length of a sampling block.The recorded samples contain the waveforms as array of amplitudes.The described data are read out scan line by scan line using the RiWaveLIB.To calculate the 3D coordinates for each sample (Section 3.3), the range r i from the sample to the point of pulse emission (origin o) is required.It is estimated as follows: The reference time (t ref ) to which the origin and direction vectors refer is derived by peak detection from the emitted waveform.For this purpose, a Gaussian distribution is fitted to the peak samples.

Multiple-time-around (MTA) processing
For data acquisition the multiple-time-around (MTA) functionality of the Riegl LMS Q680i was applied.This method allows larger measurement ranges with a high pulse repetition rate because new laser pulses are emitted before the echo detection of the previous pulses is finished.This leads to ambiguities in the assignment between emitted and reflected signal or, more precisely, the emitted and reflected signals stored in one Pulse record (Fig. 2) are not associated.
The data were acquired with MTA zone 2 what means that the echo signal of one laser pulse should be recorded after the next emitted laser pulse.Due to varying constellations of emitted and corresponding echo signals (Fig. 3), this simple assignment approach was not applicable.Rather a semi-automatic method was implemented which realises the correct assignment by using a preset range gate.The limits of the range gate r min and r max can be approximated accordingly to the flight height above ground and the height of reflecting objects.To avoid ambiguities, the range gates of successive laser pulses must not overlap.This can be achieved by setting the width of the range gate smaller than the distance a laser pulse can cover between two pulse emissions.Hence, the allowed range gate width depends on the pulse repetition rate of data acquisition.Considering the pulse repetition rate of 400 kHz and the flight height, a correct assignment could be obtained with a range gate from 500 m to 800 m.For optimisation purposes, the range gate limits are converted internally to time limits Δt min and Δt max rearranging Eq. 1.
During the data reading process, the reference time of each emitted signals is compared to the start time of each following sampling block of reflected signals (t sosbl ) while the time difference is smaller than the time gate maximum (Δt max ).Each sampling block whose t sosbl is within the time gate is assigned to the respective reference signal.

Geo-referencing
Due to occlusions which commonly occur during data acquisition of forest areas, a combined analysis of different flight strips or rather of laser pulses from different scan angles is required to achieve a detailed voxel representation of forest stands.This requires an accurate geo-referencing of the laser scanner data which will be described in the following.
The geo-referencing comprises the transformation of the origin vector (o SOCS ) as well as the rotation of the direction vector (d SOCS ) of pulse emission to a global coordinate system (Eq.2).Afterwards, the global coordinates of each waveform sample are calculated by polar point determination using the transformed pulse direction and the calculated sample ranges (Eq.1).The following right-hand coordinate systems have to be considered: 1. Scanner own coordinate system (SOCS) 2. Coordinate system of inertial measurement unit (IMU) 3. Platform coordinate system (BODY) 4. Navigation coordinate system (NED) 5. Cartesian WGS84 (ECEF) The transformations between the first three coordinate systems are platform intern and require precise information about the boresight alignment ( ).
Whereas the transformations to the navigation ( ) and to the global Earth-centred Earth-fixed (ECEF) coordinate system ( ) require access to time-solved trajectory data.

Transformation from SOCS to BODY system:
The transformation from the scanner coordinate system to the BODY coordinate system can be described with: Considering the orientation and position of each coordinate system (Fig. 4), the rotation from the scanner to the body system ( ) consists of a rotation around the y-axis and around the z-axis (Eq.4).

= (-90) • (90)
(4) After correcting alignment deviations between the coordinate systems with rotation matrix R CAL (Eq.5), the lever arm from scanner origin (o SOCS ) to IMU origin (o IMU ) is considered with the translation vector .The correction angles (Φ C ,Θ C ,Ψ C ) as well as the lever arm are related to the BODY coordinate system.Considering the navigation angles, the BODY system is rotated to the navigation coordinate system (Eq.7) which is located in the current position (λ, φ) with north-east-down (NED) orientation (Fig. 5 green).
The rotation matrix , which coincides the axis of Navigation coordinate system with the ECEF axis (Eq.8), includes the current GPS position of the aircraft (λ, φ) and can be deduced from Figure 5.The translation vector is obtained by converting the ellipsoidal WGS84 coordinates to Cartesian coordinates considering the parameters of the reference ellipsoid.

=
+ ( ) Figure 5. Relation between navigation (green) and ECEF coordinate system (blue) The corresponding emitted and geo-referenced reflected waveforms are stored in a suitable binary data format which provides the basis for further data processing.

Structure of voxel space
To realise voxel based data analysis, the geo-referenced fullwaveform laser scanning data are transferred to a 3D voxel representation.Instead of inserting discrete points extracted by Gaussian decomposition, the voxel attributes are obtained directly from the digitised waveforms (Fig. 6 left).The voxel space generation is focused on the following requirements: -memory efficient -updateable -occupied and free voxel storable To achieve memory efficiency, the voxel space is based on an octree structure which stores only relevant voxel.These can be either occupied voxel containing waveform samples or empty voxel only traversed by a laser pulse (free voxel).An octree is a hierarchical tree structure which is often used for the division of large spatial data sets.Starting at the root node which involves the whole data set, each node which contains waveform data is divided into eight sub-nodes (octants) until the bottom level (maxLevel) of the octree or rather the desired resolution of the voxel space is reached (Fig. 6 right).There are two node classes to be defined: the leaf nodes on the bottom level which represent the final voxel and hold the voxel attributes and the parent nodes on the previous levels which hold a pointer to their eight child nodes.The size of the root node (size octree ) is adapted based on the extent of the full-waveform data (size data ) and the specified resolution of the voxel space (Eq.9).To allow the consideration of laser scanner characteristics (e.g.footprint and sampling interval), the horizontal (vs XY ) and vertical voxel size (vs Z ) can be defined independently.
Depending on the chosen voxel size or if multiple flight strips are inserted into the voxel space, a voxel can be hit several times.To allow an updating of the voxel attributes, each voxel stores a list of values required to obtain the final voxel attribute (Section 4.3).In addition to generate new voxel, specific voxel can be deleted which is essential for the filtering process.

Integration of waveform data
After determining the bounding box of the root node and the number of needed levels, the octree can be build by intersection with the waveform data.In this process, the data are interpreted as ray (x r ) with amplitudes at specific ranges (r).Thus, for each laser pulse emission, a ray can be defined considering the georeferenced origin and unit length direction vector (Eq.10).

= + • (10)
On this basis, the implemented ray tracing approach was developed combining the parametric algorithm presented by Revelles et al. (2000) and the voxel traversal method proposed by Amanatides & Woo (1987).Both methods are well suited for this application because all computations and comparisons are related to the ray parameter r.The first approach is applied for checking the intersection between the ray and a current node including the detection of the entry and exit plane and detection of the first crossed sub-node.Beyond that an adaption concerning negative components of the direction vector was necessary.In case of dividing a node, the traversal through the eight sub-nodes was realised interpreting the sub-nodes as 2 x 2 x 2 grid (Fig. 7).The coordinates of the eight sub-nodes are defined as binary vectors and are the basis for the traversal method of Amanatides & Woo (1987).Because the implemented ray tracing approach is a top-down method, the intersection of the waveform data ray with the voxel space begins at the root level.Rearranging Eq. 10, the ray parameters and which correspond to the intersection of the ray with the root boundaries (Fig. 8) are determined where index i defines the coordinate direction (0 ≙ x, 1 ≙ y, 2 ≙ z-direction).Using Eq. 11 the range interval within the node boundaries is obtained.Each waveform sample k for which r min < r k ≤ r max is within the node.
If r min < r max the ray intersects the root node and it is added as parent node in the octree structure.Additionally, every detected parent node is appended to a node list which allows an iterative setup of the octree.In the division process, at first the ray parameters , which correspond to the boundaries dividing the parent node one half each (Fig. 8), are calculated (Eq.12).These values are used to establish the minimum and maximum ray parameter of the subnodes.
The planes, where the ray entries and exits a node, are defined by the index of the maximum or minimum ray parameters (Tab.1).

Maximum Entry plane Minimum
Exit plane YZ YZ XZ XZ XY XY Table 1.Determination of entry and exit planes of a node Considering the relation of and r min , the coordinates of the first crossed sub-node are obtained.The detected exit plane of the sub-node defines the coordinate direction of the next subnode along the path of the ray.The coordinates of the next subnode are obtained by adding a step to the corresponding node coordinate.The step is 1 if the respective component of the ray direction vector is positive and -1 if it is negative.If the ray exits the parent node, the next parent node from the node list is processed in the same way.Each intersected node is analysed whether it is a parent node, an occupied leaf node or a free leaf node.If the node list is depleted, all voxel crossed by the ray are created as leaf nodes in the octree structure.The final voxel space is generated if every waveform ray was intersected with the octree.

Determination of voxel attributes
After integrating all waveforms into the voxel space, redundant voxel entries especially occurring in overlapping areas of different flight strips have to be taken into account.Hence, voxel attributes shall be determined from the waveform amplitudes.To consider attenuation effects caused by partial reflections in higher parts of the canopy, the amplitudes of the raw waveforms were corrected applying the integral method published by Richter et al. (2014) based on the reference value for leaf-off data proposed by Richter et al. (2015).This allows the calculation of radiometrically corrected voxel attributes for forestry applications.
To emphasise pulses which are closer to the nadir, two scan angle dependent methods for the determination of the voxel attributes (VA) were implemented: (1) use the amplitude corresponding to the smallest scan angle and (2) calculate a scan angle dependent weighted average value (Eq.13).To realise this, the attribute lists of the occupied voxel were filled with amplitudes (a) and corresponding scan angles (ϑ) of the contained waveform samples.The calculation of a weight (w) is based on the ratio between scan angle of the laser pulse and the maximum scan angle of data acquisition (ϑ max ).
= number of voxel entries ⁄ (13) Both methods were compared to voxel attributes resulting from the maximum amplitude of a voxel.In particular, differences in the inner parts of the canopy (Fig. 9 a) and along the tree stem (Fig. 9 b) are obvious for both methods.The large differences between the methods (Tab.3) reveal that smaller scan angles do not necessarily correlate with a smaller footprint and a higher backscattering.In fact, the orientation of the reflecting object surface and the resulting incidence angle have a great influence.
angle minimum 0] weighted mean value 1] Table 3. Ranges of differences between scan angle dependent methods and maximum amplitude In addition, due to integrating unfiltered waveforms in the voxel space, the noise influences the averaging and leads to decreased attributes.Further analysis is required to evaluate the suitability of a scan angle dependent determination of voxel attributes.To ensure the comparability of the following filtering methods, the maximum amplitude is initially preferred as voxel attribute.

Method
The resulting voxel space contains a lot of voxel which do not correspond to object information.To reduce the voxel number and to extract only relevant object information, a filtering process is followed.Instead of filtering the waveforms, a 3D filtering of the voxel space is performed which allows the consideration of neighbouring information.Concerning the vertical orientation of trees, a Gaussian filter kernel is applied column-wise to the voxel space to calculate temporally new attributes for each voxel.Based on these attributes, each column is analysed applying an adapted thresholding.Beginning at the top, the attributes are successively compared to a threshold.Voxel with a value smaller than the threshold are marked as non-object and above the threshold as object voxel.This strict proceeding yields a cutting of the waveforms and a loss of information.To keep the entire peak form, non-object voxel adjacent to object voxel are remarked as object while their attributes are descending.Finally, all non-object voxel are deleted from the voxel space.

Determination of threshold
The used threshold is obtained by a histogram analysis of all voxel attributes.Lefsky et al. 2005 proposed a histogram-based approach to determine waveform dependent thresholds.This is working well in case of a small signal-to-noise ratio of the waveform.Otherwise the threshold determination is more uncertain and tends to larger thresholds.
Considering the whole voxel space, the voxel attributes show a great signal-to-noise-ratio and a clear differentiation of nonobject and object voxel in the histogram.Hence, this approach is well suited to obtain a voxel space dependent threshold.At first the histogram of the voxel attributes is generated.
Assuming that the most non-object voxel are represented by the maximum peak, a Gaussian function is fitted to the respective bins.The threshold results from the resulting mean plus three times the standard deviation.

Geo-referencing
The assessment of the geo-referencing was realised comparing point clouds at reference surfaces.Indicating that the detection of suitable reference surfaces is difficult in forest areas, 12 building roofs and one road were used to calculate cloud-tomesh distances.For this purpose, point clouds were generated from the geo-referenced waveform data by fitting a Gaussian distribution to detected peaks.To assess a relative accuracy, the flight strips were compared in overlapping areas, whereat strip five acts as referee and as basis for the Delaunay triangulation.Considering the mean values of the signed distances, the comparison of the flight strips yield differences of a few centimetres.Though it shows that the data reading and transformation was implemented correctly, it is obvious that there are still errors (e.g. of boresight alignment) which have to be eliminated by strip adjustment.In particular the flight strips one and two are shifted to strip five nearly in the same direction.This is also obvious from Figure 10 which visualises exemplarily the signed cloud-to-mesh distances of one building roof for flight strips one to four.

Voxel space generation
According to the results of geo-referencing, a voxel space of a spruce was generated by integrating the waveforms of flight strips three, four and five.As flight strip five provides the nadir direction (about -3.7°), the other two strips yield scan angles of about -17.5° and 28.8°.Considering the sampling interval and the footprint of the laser scanning data, the voxel size was set to 30 x 30 x 15 cm³.The extent of the selected data set yield 12 x 25 x 43 m³ (size data ).Hence, to achieve the predefined voxel space resolution, nine divisions (maxLevel) are required (Eq.9).The initial extent of the voxel space (size octree ) results to 153.6 x 153.6 x 76.8 m³ (Eq.9).The selected data set comprises 1726 pulses which were intersected with the voxel space whereat only occupied voxel were registered.Figure 11 shows the generated filtered voxel space as well as the corresponding point cloud obtained by Gaussian decomposition for comparison reasons.The generated voxel space contains 88674 voxel (Tab.5) whereat voxel with one entry (Fig. 12 right green) provide the majority (68 %).

Voxel space filtering
Figure 12 shows left hand side the relevant part of the histogram derived from the attributes of the unfiltered voxel space overlaid with the fitted Gaussian distribution.A resulting threshold of 8.002 was determined.The 3D filtering was performed applying a 3 x 3 Gaussian kernel.Hence, for each voxel the surrounding area of 90 x 90 x 45 cm³ was considered in the filtering process.In addition to the 3D filtering, a waveform based filtering proposed by Persson et al. (2005) was performed for comparison reasons.By this, each waveform was filtered applying a waveform dependent threshold which was defined with three times the standard deviation of noise.Afterwards, the filtered waveforms were integrated into a voxel space with the same parameters.Table 5 opposes the results of the filtering methods for some parameters.It is obvious that the 3D filtering eliminated more voxel (57%) than the waveform filtering.The unchanged ranges of the voxel attributes prove that the adapted thresholding is working well and the peak samples below the threshold are kept.The histograms of the redundancies depict that most deleted voxel contained only one entry (Fig. 12  To achieve a more detailed analysis of the filtering methods, Figures 13 and 14 visualise comparatively the results with the help of a vertical slice and a profile of a column extracted from the unfiltered and filtered voxel spaces.It is obvious that the 3D filtering performs better than the waveform filtering because more noise voxel could be eliminated whereas the object parts especially the canopy voxel are kept.Nevertheless, there are still remaining noise voxel below the ground which have to be eliminated.For instance, this could be obtained by detecting the digital terrain model from the voxel space and neglecting all voxel below it.

CONCLUSION
The paper specified the entire procedure from accessing the raw full-waveform laser scanning data to the successful generation of an attenuation-corrected and filtered voxel based representation of the waveform data.To achieve a correct data reading, a semi-automatic MTA processing algorithm was described which is applicable to different MTA zones.
Currently, only laser scanner data acquired with a constant MTA zone are analysable.The implemented ray tracing approach is a valuable option for the efficient integration of ALS waveform data into a voxel space and allows additionally the analysis of the pulse history.The major advantage of using an octree as data structure is the memory efficiency.The successfully generated voxel based representation of a spruce results from the combination of three flight strips.Considering redundant voxel entries occurring in the overlapping areas of the flight strips, the comparison of the maximum amplitude with scan angle dependent voxel attributes shows decreasing values in parts of the inner canopy and the tree stem.Hence, further analysis is required to evaluate the suitability of these methods for forestry applications.
The results show that the developed 3D filtering approach of the voxel space performs better than the commonly used waveform filtering.The specified method for voxel space generation can easily be adapted to the integration of the differential backscatter cross section which will be realised in the next step.
For further enhancement of the volumetric representation of forest stands, future work will deal with the consideration of the laser beam divergence in the voxel space generation.

Figure 1 .
Figure 1.Generalised structure of SDF file format

Figure 2 .
Figure 2. Structure of Pulse record

Figure 3 .
Figure 3. Principle of semi-automatic MTA processing

Figure 4 .
Figure 4. Orientation of platform based coordinate systems (left to right): scanner, inertial measurement unit and platform coordinate system 3.3.2Transformation from BODY to ECEF: To convert the coordinates from the BODY system to the global coordinate system, two remaining transformations are necessary (Bäumker & Heimes, 2001).The required position and orientation of the aircraft has to be interpolated in the timesolved trajectory data.This binary POF file contains the ellipsoidal latitude, longitude and height (λ, φ, h) from GPS measurements (WGS84) referred to the IMU origin and the navigation angles roll, pitch and yaw (Φ, Θ, Ψ) from IMU measurements with 3.9 msec resolution.The correct instant of time for interpolation is achieved by converting the internal reference time (t ref ) to GPS time (Eq.6).= + − (6)

Figure 6 .
Figure 6.Principle of ray tracing (left) and corresponding generalised octree structure (right)

Figure 7 .
Figure 7. Numbering and coordinates of sub-nodes

Figure 9 .
Figure 9. Two vertical slices of the voxel space coloured by absolute differences between amplitude of minimum scan angle and maximum amplitude

Figure 10 .
Figure 10.Meshed building roof (left) and corresponding histogram of signed cloud-to-mesh distances (right) for flight strip one, two, three and four (top to bottom) The flight strips three, four and five show the best match.The other two flight strips are shifted about 4 to 5 cm downwards approximately along the right roof area.The roof points of strip three show a small shift which is not obvious from the road points.Additionally, an absolute accuracy is estimated comparing the geo-referenced full-waveform point cloud of strip five to a conventional airborne laser scanning data set acquired 2006 with 1.5 pts/m² offered in UTM coordinates with DHHN92 heights.Thus, a coordinate transformation of the full-waveform point clouds was required.Calculating the signed distances of the conventional ALS data to the meshed point cloud of strip five yields a mean value of -0.05 m.Hence, the full-waveform data are located above the conventional ALS data.The large standard deviation of 0.057 m is caused on the one hand on the dispersion of the ALS points and on the other hand on a remaining coincidence error of both point clouds.This can be effected by the transformation of the height coordinate systems as well as by remaining boresight alignment errors in the georeferencing process.Nevertheless, the dimension of the discrepancy between the data sets confirms the correct proceeding of the applied geo-referencing.

Figure 11 .
Figure 11.Visualisation of a spruce coloured by attenuationcorrected amplitudes: Point based (a) and voxel based representation (b) of waveform data

Figure 12 .
Figure 12.Histograms of voxel attributes (left) and voxel entries (right) of unfiltered and filtered voxel space

Figure 13 .
Figure 13.Results of waveform filtering: vertical slice of unfiltered and filtered voxel space overlaid (a), vertical slice of filtered voxel space coloured by attributes (b), profiles of a column (c)

Table 2 .
Table2summarises the properties and further procedure for each node type.If a node does not exist in the octree structure, it is added as new subnode to the respective parent node considering the node type and the node coordinates.Properties and procedure for the different node types Table 4 summarises the statistical parameter of the signed distances.

Table 5 .
left).Parameters of voxel space before and after filtering