9.4.3 LIDAR Interpolation

In general, there are two principal ways to prepare LIDAR data for filtering techniques: gridding or working on the original point cloud. On the one hand, some researchers suggest the use of the original data for accuracy reasons, since gridding and interpolation may omit original data and induce additional non-existent data. On the other hand, gridding and interpolation allows the use and further advancement of well-known image processing techniques [47], such as filtering in the spatial and frequency domain. For applications requiring data fusion, different bands have to be co-registered to obtain a consistent format. Typical spatial resolutions for gridding are in the order of 0.25–2 m [104]. Commonly used interpolation techniques include nearest neighborhood and bilinear interpolation [34, 111, 146, 167, 181].

Bater et al. [20] compared seven interpolation methods applied to LIDAR data at different resolutions and measured the interpolation error as RMSE. Conclusions were that, although the interpolation error is independent of the interpolation method, it is a function of the resolution, where 0.5 m increase in spatial resolution yields an RMSE increase of 1 cm on average. Therefore, the authors suggested employing a natural neighbor algorithm with minimal interpolation error and computational costs. Natural or nearest neighbor interpolation suggests that missing data can be synthesized by the mean (or median) of adjacent data points based on the assumption that they are similar.

9.4.4 LIDAR Filtering

The main interest in airborne LIDAR data is to generate accurate elevation maps and to identify objects within the point cloud. The main products include Digital Surface Models (DSM), Digital Terrain Models (DTM) and normalized DSMs (nDSM) [143]. A DSM includes all sampled top surfaces of objects and ground [14, 136, 138], including tree canopies, roof tops, chimneys or above ground power transmission lines, as illustrated in Fig. 9.12. A DSM can also be directly produced from the LIDAR point cloud. Since many applications require accurate ground data, one major goal in LIDAR filtering is to separate object and ground points. If object points are of interest, an nDSM can then be generated very conveniently by subtracting the DTM from the DSM [172]. Again, an interpolation step is required to create a DEM from the filtered data (to fill those patches where objects were with appropriate data).

LIDAR filtering is the bottle-neck in the processing chain from acquisition to real applications. Friess [50] reported that the ratio of post-processing to acquisition time for correcting the points’ geometry is up to 14:1 and DEM generation demands an even higher ratio of 20:1. Filtering algorithms constitute the major contribution to a final product derived from LIDAR [167], taking three times longer than the acquisition. Flood [43] reported that at least 60 % of the production time to obtain DEMs is allocated for filtering and other post-processing steps. Although LIDAR

Fig. 9.12 From top to bottom: Original scene, DSM, DTM, nDSM

systems are considered to be mature, there is much work to de done to develop filtering algorithms that are robust, accurate and efficient [10, 13, 162]. In the following, a brief overview of the different categories of LIDAR filtering techniques is given.

Morphological filtering does two things: filters object points and interpolates the missing terrain at the same time. Based on set theory, morphological filtering consists of two basic operators, erosion and dilation [65, 66], usually applied in a sliding window. Further derived operators are cleaning, filling, bridging and the watershed algorithm. When applied to LIDAR data, the major assumption is that there is a distinct difference between slopes of objects and of ground [172]. The first step is to find the minimum LIDAR point in a window [90] where the difficulty here is to select the right window size according to shape, size and orientation of objects [8]. Small windows relative to large objects result in treating them as ground. Large windows, however, may erode terrain irreversibly. Therefore, the latest developments in

morphological filtering make use of progressive [181] increasing window sizes [90]. In general, morphological filtering is considered to be robust and works well for isolated objects. However, some drawbacks include the indifference of this filtering technique towards object and terrain details, which can be eroded irreversibly [32].

Slope-based filtering exploits slope and discontinuities between adjacent single or grouped points (e.g. buildings) since LIDAR data is vertically accurate in the range of decimeters. Vosselman [155] employed a predetermined cut-off plane to measure and classify slopes of lines between two points with parameters to adapt the plane to highly sloped terrain [142]. Slope-based filtering is non-trivial since extreme slopes and flat terrain can cause failure of the algorithm. Therefore, the Laplacian of Gaussian (LoG) is often used to differentiate objects from the ground [55]. A grid based DTM generation approach was developed by Wack et al. [163, 164] who downsampled the data coarsely first to a 9 m spacing. By doing so, most of the object points were discarded. Roggero [128] derived a DTM by using weighted linear regression on original, irregular LIDAR data. The authors then employed a local inverted cone adapted to the terrain slope in order to estimate ground points, regarding the maximum height difference between two points. Objects were then detected using a LoG. In general, slope-based algorithms pick up precisely edges, lines, corners and other discontinuities. If spatial masks are used, gridding is required. Slope-based filtering algorithms fail if their major assumption steep slopes are only caused by objects is violated; for example, in mountainous terrain with discontinuities such as trenches, manhole covers for canalisation, holes, caves or outliers.

Geometry-based filtering is a popular method to describe man-made and natural objects in LIDAR data, assuming that object properties such as shape, length, width, height, position and orientation are known. From an nDSM obtained by means of morphological filtering, Weidner and Förstner [172] extracted buildings from a parametric and prismatic building model [169] based on the Minimum Description Length (MDL) principle [170]. First, a bounding box for each building is estimated using a defined threshold with respect to known building heights. Second, the segmentation is refined based on a height threshold calculated from the median of 10 % lowest and 10 % highest points of the preselected data. Sparse Triangular Irregular Network (TIN) densification for DTM generation is also very popular [145] since neither gridding and interpolation nor dense point clouds are required [96]. Axelsson [10] classified vegetation, buildings and power transmission lines in original LIDAR data, based on their surface in a TIN using MDL criteria. The underlying assumption was that objects consist of planar faces and therefore, neighbouring TIN facets had a similar orientation with a second derivative of zero whereas the second derivatives of breaklines and vegetation points were non-zero. One of the biggest advantages of geometry-based filtering is that structures in LIDAR data can be described directly from point clouds. Buildings can be recovered very efficiently and initial gridding and interpolating the data is not necessary. However, geometry-based filtering requires many predefined parameters and becomes more difficult when the level of object complexity is high.

Curvature-based filters have been developed to detect classes of curved areas (i.e. convex, concave and plane) within the point cloud [6, 153]. Vosselman [154]

assumed that buildings consist of planar faces which can be recognized by applying the Hough transform in LIDAR point clouds. The data had been segmented with existing 2D ground plans [156, 157]. The spatial structure chosen for the geometric description of the planar faces were 2D Delaunay triangles [154]. If a connected component was greater than a specified threshold, a planar face was found. The buildings’ outlines were estimated using the planar faces projected into the 2D space [154]. Working on non-gridded data, Rottensteiner and Briese [131] detected non-ground building regions in the first stage. Then, roofs were detected using surface normals of a DSM [134]. The biggest advantage of curvature-based LIDAR filters is that they are directly applicable to the original data to detect structures in a LIDAR point cloud. Drawbacks, however, are that only man-made objects that have planar, convex or concave faces can be detected with pre-defined thresholds and that vegetation is difficult to recognize due to its random character.

Linear prediction can be used for DTM generation and gross error removal [26], as was developed by Kraus and Pfeifer [93]. The supervized algorithm contained a combination of filtering (object point classification) and interpolation of the ground. First, a rough approximation of a DTM was calculated using linear prediction employing overlapping patches. The difference of this DTM and the original LIDAR point cloud was estimated to obtain residuals as the basis for assigning weights to LIDAR points from a parametric weighting function. Points above a certain threshold were classified as object points. The whole algorithm was iteratively executed until either a stable situation or a predefined number of iterations was reached [93]. Based on linear prediction, Kobler et al. [91] addressed the challenge of filtering LIDAR data in steep wooded terrain with a repetitive interpolation algorithm by setting empirically estimated parameters, thresholds and a buffer zone [94]. A major advantage of linear prediction models is their applicability to different terrain types and that DTMs can be derived from sloped terrain. The compromise to this solution however is that prior knowledge of the terrain for the weighting factors and thresholds are required. Those parameters can be adjusted depending on the application and terrain. Also, some approaches lack a clear termination criterion, as either a stable situation has to be reached or the number of iterations has to be specified in advance.

Coarse-to-fine filters and multi-resolution filters benefit from both detailed and coarse views of the LIDAR data [93]. Pfeifer et al. [120] presented a series of LIDAR post-processing steps aimed at flood modeling: data calibration, strip adjustment, ‘robust linear prediction’ and terrain structure recovery (i.e. breakline modeling). The new element contributing to the robustness of the filtering approach was the use of hierarchical pyramids. DTMs of different resolutions were compared to each other, however, still requiring some thresholds [120]. References [139–141] presented a semi-automatic, multi-resolution approach for filtering LIDAR data based on the Hermite transform. Gaussian pyramids were employed to transform the data in the first step. An assumption was made that any change from ground and off-terrain is a linear combination of shifted, rotated and scaled error functions, similar to a kernel function in Wavelet analysis. A multi-resolution algorithm for DTM generation from LIDAR data was proposed by Vu et al. [159–162] who compared

successive median-filtered resolutions of gridded LIDAR data to detect boundaries. The final segmentation of the LIDAR data was achieved using both the boundaries and the actual height as features. A Wavelet approach to separate ground and object points on gridded LIDAR data was proposed by Vu and Tokunaga [158]. The authors applied K-means on height to segment buildings, motorways, boundaries and trees. The advantage of multi-resolution filtering algorithms is the separation of low and high frequencies (i.e. approximations and details) of a LIDAR scene. However, it is not always clear what resolution to choose since energy (i.e. information contents) becomes smaller the further the signal is decomposed [18]. When using approximations (i.e. discarding the high frequency components) a loss of information occurs. Also, higher computational costs and memory requirements have to be considered, when using multi-resolution filters.

Knowledge-based filters make use of predefined models of different height and shape [22]. Haala and Brenner [60] presented a complete solution of a realistic 3D city model from LIDAR derived DSMs. To get the buildings’ ground planes, the authors segmented the DSM using a 2D GIS map whose incompleteness and inaccuracy was complemented with a cadastral map. Their algorithm was based on region growing using straight lines of pixels [59] and histogram analysis of surface vector normals which yielded planar surfaces [62]. Having obtained the ground planes, the actual buildings were modeled by fitting a limited number of predefined 3D building primitives into the DSM. Location and rough type of vegetation were estimated with three bands of aerial Color Infra-Red (CIR) images (i.e. NIR, red and green) and the geometric information from the nDSM [63]. A hierarchical rule-based filtering solution was presented by Nardinocchi et al. [114]. LIDAR data was gridded and three classes—terrain, building and vegetation—were estimated by exploiting their geometric (height differences) and topological (spatial distribution) properties and their relationships among each other using region growing and local slope analysis. Employing hierarchical rules and fitting in 3D geometric primitives is very efficient for identifying objects in LIDAR data [114], provided that urban area is presented. However, this approach involves prior knowledge of the buildings and other manmade structures and becomes extremely difficult for complex objects. Furthermore, it is almost impossible to model vegetation geometrically. Hence, additional information such as CIR or Normalized Difference Vegetation Index (NDVI) has to be integrated.

Data fusion for land cover classification exploits complementary properties of both LIDAR and photogrammetry. It requires registration and orientation of gridded and interpolated LIDAR data with all spectral bands [136, 137]. Popular spectral data are NIR, red and green from CIR imagery [61] to be combined with height information from LIDAR data. Processing co-registered remotely sensed data involves solving a multivariate statistical problem. Typical classifiers are used such as ML, distance classifiers, Support Vector Machine (SVM), Principal Component Analysis (PCA), Independent Component Analysis (ICA) and Artificial Neural Networks (ANN). CIR imagery combined with LIDAR data can be used for detecting sealed and non-sealed surfaces for water waste management and council taxing [138] or for tree crown volume estimation. Further applications with fused LIDAR data and

aerial images are building reconstruction [131], building outline extraction with fused NDVIs from IKONOS multispectral panchromatic imagery [144], roof segmentation [23, 24], building detection based on Dempster-Shafer [132, 133], roofs, trees and grass classification using Gaussian Mixture Model (GMM) and Expectation Maximization (EM) [30], per-pixel minimum distance classification of streets, grass, trees, buildings and shadows [61], including empirical ground truth collection [118]. For data fusion applications, it is popular to use open source software such as the Geographic Resources Analysis Support System (GRASS) [27] or commercial software packages such as eCognition, Arc/Info, TerraSolid and TerraModeler [76, 152, 171]. In general, fusing LIDAR and additional space and airborne remotely sensed data has great potential for improving accuracy for land cover classification as it combines advantages of complementary bands. The gain, however, has a limit due to the curse of dimensionality in that an extensive adding of further bands may result in attenuation of accuracy [83]. Furthermore, it is non-trivial to find contemporary 2D GIS maps, imagery and ground truth. Moreover, the data has to be co-registered, and if necessary downsampled, gridded and interpolated due to different resolution and orientation. This process, however, omits data and adds additional non-existent information. The use of commercial software packages involves experimental tuning of the settings and most algorithms are hidden to the user. Furthermore, expensive licenses cannot always be purchased due to funding limitations.

Statistical classification filters are a means to segment and classify objects and terrain in LIDAR data in an unsupervized manner [168]. Cobby et al. [34] generated a DTM for flood simulation from LIDAR data recorded by the Environment Agency, UK. The authors segmented the rural area close to the River Severn and classified the objects into three vegetation height classes: short (crops and grasses), intermediate (hedges and shrubs) and tall (trees). Using a 5 × 5 window, the semiautomatic segmenter employed the standard deviation of local height as a feature because tall vegetation was assumed to constitute the highest objects [35]. For simplicity, the authors declared the limited number of houses in this area as high vegetation. Since the authors were interested in vegetation height, they separated the slightly hilly terrain from the actual object using detrending [39]. By subtracting the bilinear interpolated surface from the LIDAR data, an nDSM for a hydraulic model was obtained, as developed by [33]. Skewness balancing [19] and its adaptation to hilly terrain [17] are an alternative to object-based filtering. Yao et al. [175] have adapted skewness balancing to address the challenging task in detecting cars in LIDAR data. Bao et al. [15] further improved the algorithm by incorporating the measure of kurtosis, where the authors exploited the different changes of both statistical measures with respect to vegetation and ground. Advantages of statistical filtering algorithms are that they can work unsupervized directly on the original, non-gridded data. They fail, however, if the statistical model boundaries are insufficiently described.

A comparison of LIDAR filtering methods is illustrated in Table 9.5.