Correction of Erroneous Model Key Points Extracted from Segmented Laser Scanner Data and Accuracy Evaluation

Abstract

Point cloud data (i.e., LiDAR; Light Detection and Ranging) collected by Airborne Laser Scanner (ALS) system is one of the major sources for surface reconstruction including DEM generation, topographic mapping and object modeling. Recently, demand and requirement of the accurate and realistic Digital Building Model (DBM) increase for geospatial platforms and spatial data infrastructure. The main issues in the object modeling such as building and city modeling are efficiency of the methodology and quality of the final products. Efficiency and quality are associated with automation and accuracy, respectively. However, these two factors are often opposite each other. This paper aims to introduce correction scheme of incorrectly determined Model Key Points (MKPs) regardless of the segmentation method. Planimetric and height locations of the MKPs were refined by surface patch fitting based on the Least-Squares Solution (LESS). The proposed methods were applied to the synthetic and real LiDAR data. Finally, the results were analyzed by comparing adjusted MKPs with the true building model data.

1. Introduction

Point cloud data obtained from Airborne Laser Scanner (ALS) system is composed of dense 3D coordinates. The key of the data process is how to organize irregularly distributed data to get meaningful information to reconstruct various terrain features. Sequence of the procedures including “filtering”, “feature extraction”, “classification”, “segmentation”, “grouping”, and “modeling” for natural and artificial objects are involved. However, a standard procedure and methods required in each step are not available yet. Instead, many of the methods are based on ad hoc approaches that sometimes yield inconsistent results due to data dependency problem. Human being could easily identify features and extract information by visualizing such point cloud data. Therefore, automation of the Light Detection and Ranging (LiDAR) data processing belongs to high-level process such as visual perception, pattern recognition, object identification, computer vision, etc.

The central issue in implementation of the processing methodology is to make algorithms resembling human capability of object recognition. Another important issue is evaluation of the results. The procedure needs to take into account trade-off between automation and quality of the results. Especially, quality assessment is based on both positional accuracy and Level-of-Detail (LoD) of the objects. Basically, LiDAR data is list of the numerical coordinates (X, Y, Z) of the points on the terrain surfaces. Therefore, data and model should be differentiated. Digital elevation data is simply collection of points having 3D coordinates, while quite sophisticated processes are required to generate Digital Elevation Model (DEM). Object modeling means reconstructing various kinds of the individual feature. Recently, Digital Building Model (DBM) and Building Information Modeling (BIM) are one of the rapidly growing applications in geospatial information for 3D urban and indoor modeling and navigation (Dongzhen et al., 2009; Kolbe, 2009).

Identification of geometric characteristics of objects and feature collection are significant tasks in modeling. In conventional photogrammetry, 3D digitizing Model Key Features (MKFs; e.g., corner points and object boundaries) is performed on the oriented stereo images. Photogrammetric method utilizes visual information based on the brightness values of the images. LiDAR data processing, on the other hand, makes use of the height information. An obvious advantage of using LiDAR data is to derive more explicit geometric information of the objects than using images. The method introduced in this paper starts with the MKPs extracted from segmented patches. Many studies on LiDAR data segmentation and modeling have been performed during the last years including Habib et al. (2006), Kim et al. (2009), Sampath and Shan (2008), Vosselman and Dijkman (2001), Vosselman and Klein (2010), Yoo et al. (2012) and so on.

How to segment and to extract MKPs are not the central issue of this study. However, the MKPs sometimes are not determined correctly due to the imperfectness of segmentation and/or detecting the model key features. Correction of the erroneous MKPs is core of this study.

Procedure of the proposed method is as follows: Surface fitting segmented data of each building object Comparison of MKP coordinates with fitting results Comparison of MKP coordinates with original LiDAR data Comparison of MKP coordinates with true building model data Correction of errors (or discrepancies) by analyzing above comparisons Evaluation of results and quality assessment

The proposed scheme was applied to the simulated airborne LiDAR data and real data with various roof types such as flat, gable, hip, pyramid, dome, half-cylinder, etc. In order to correct coordinates of the MKPs, surface fitting was performed by least-squares with all of the segmented data, and new heights of the MKPs were computed using the coefficients of the surface equations. In addition, instead of using all segmented data, part of the data (i.e., central region of the segmented data excluding around the segment boundaries) were used to derive surface equations since most probably some errors might occur along the boundaries. New heights were also computed in this case. If the height value differences between original MPKs and computed MPKs were larger than the tolerance, the original height values were replaced by computed ones. Basically, the discrepancy at each MKP was compare with the true data for the analysis and evaluation.

Overall concept of the segmentation and problems are described in Section 2. Afterwards, determination of the MPKs from the segmented surface patches, and correction method of the MPKs are introduced in Section 3 and 4, respectively. Results from the experiments and analysis are shown in Section 5. Conclusions and recommendations for the future work are presented in the final Section.

 

2. Surface Patch Segmentation

2.1 Overview of Segmentation

Segmentation is a process of identification, partitioning and grouping data into meaningful and homogeneous regions or structures. The ultimate goal of photogrammetry is to create model world (e.g., city/building models, digital topographic maps, etc.) from various spatial data including aerial/satellite imagery or LiDAR data. Sometimes, two or more different sources of the data are integrated. The model world represents abstract of the real world and contains representative information. Therefore, standardized feature collection rules such as generalization and regularization of the reconstructed world have to be applied to generation of the model world (Schenk, 1999). In automatic object modeling using LiDAR data, surface patch segmentation is one of the essential steps because elements of the features (e.g., MKPs) to depict objects could be extracted from the segmented surfaces. However, results of the segmentation are subject to characteristics of the data and processing method.

Fig. 1 shows the typical workflow of airborne LiDAR data processing for object modeling. Since this paper focuses on the evaluation and correction of the MKPs extracted from segmented patches, the preceding steps such as preprocessing, filtering and classification are not the main issues in this study. Each step in the procedure affects quality of the modeling eventually.

Fig. 1.Typical workflow of object modeling using airborne LiDAR data

One of the most common segmentation methods of the LiDAR data is those that group points that fit to the same surface such as plane, smooth or curved surfaces (e.g., sphere, cylinder, etc.). Shapes in the point cloud could be recognized based on the segmentation. In most case, planar segments are required for modeling buildings. Various methods of segmenting point cloud data are available such as 3D Hough transform, Random Sample Consensus (RANSAC), surface growing, and scan line segmentation (Vosselman and Klein, 2010; Li and Guan, 2011). Lim (2008), and Park et al. (2012) utilized 3D ψ-s curve, and chain code methods to segment airborne LiDAR data for building modeling, respectively. These two studies are based on the shape descriptors that are used for object recognition in image processing and computer vision.

The fundamental elements of the segmentation include proximity, direction, slope, continuity, collinearity, similarity, closure among the points. Some of the factors quite coincide with Gestalt laws of grouping (or organization) in highlevel visual perception (Ullman, 2000). The challenging work in segmenting LiDAR data is to group points into the same surface patch automatically without visual analysis. Since LiDAR data basically provides 3D coordinates of the points, to derive all necessary elements with appropriate condition and constraint is difficult problem. Therefore, the segmentation algorithms have to be adaptive approach on the bases of data properties such as point density, topography, morphology, and the like. Moreover, the essential properties of the segmentation scheme are to be scale, translation and rotation invariant.

2.2 Problems of Segmentation

Often the segmentation boundaries are not clearly defined due to characteristics of the data and shape of the objects. The common problems are over/under segmentation, and invading/invaded segmentation that might not be avoidable in many cases as shown Fig. 2. For examples, a roof plane is segmented into several patches, or adjacent roof planes with different slopes are segmented as one patch. Therefore, segmentation is matter of the quality assurance (Habib, 2013). Fig. 3 demonstrates typical problems of segmentation. Patch A includes data from ground and other patch data, patch B is disturbed by trees, and patch C and D have overlap regions each other.

Fig. 2.Examples of erroneous segmentation

Fig. 3.Typical problems of segmentation

There are several reasons why the segmentation results often are not satisfied: i.e., noises in the raw data, density of the points, other objects near the buildings (e.g., trees, cars, etc.) interpolation if applied, imperfection of the algorithms including preprocessing and segmentation, and so on. Even though noise removal, filtering for ground/non-ground feature separation, and building extraction are preceded before segmentation, the methods involved with these procedures may not be perfect. Therefore, it is worthwhile to find methods to provide accurate MKPs for building modeling under such situations.

2.3 Segmentation Method with 3D Chain Code

In this study, segmentation results were directly adopted from the multi-cube with 3D chain code method published by Park, et al (2012). The concept of this method is that 3D slopes of the point clouds are converted to the chain codes and grouped based on the codes. For more detail description of the shape, sub-cubes are generated by analyzing distribution of the points. Scheme of the method is depicted in Fig. 4., and brief introduction of the method is as follows: Create 3x3x3 cube with chain codes, i.e., 3 layers to compute 3D distribution of the points Generate sub-cube based on the point distribution, i.e., adaptive hierarchical approach Apply the cube for 4 directions to obtain rotation invariant chain codes Grouping of the chain codes for each direction Combine the grouping results for final segmentation

Fig. 4.Segmentation of point clouds with multi-cube chain code

 

3. Model Key Points

3.1 Role of MKPs

One of the crucial procedures involved with object modeling is feature collection. In photogrammetry, 3D digitizing of model key features such as points, lines and polygons of the objects on the stereo images are performed. In general, extracting model key features (e.g., boundaries of the patches, and corner points of the patches) is followed by segmenting surface patches for object modeling. MKPs have to be well-defined and distinctive points that represent overall characteristics (i.e., shape, size, slope, orientation, structure, etc.) of the objects. However, the model key features could not be accurate due to incompleteness of the segmentation. It is obvious that quality of the segmented patches affects accuracy of the model key features.

The frame of a building object can be formed by simply connecting between MKPs, so called, CityGML LoD-2 building model. LoD-2 is block structures combined with roof features. LoD-2 building models can nowadays be created efficiently using airborne LiDAR data. Therefore, derivation of the geometrical features of building roofs is one of the main research topics of the LiDAR based 3D building modeling (Boeters, 2013; Straub et al., 2009). Fig. 5 shows building models of the LoD-2. Fig. 5(a) is simulated buildings used in this paper. The LoD of the objects is related to the cartographic generalization as shown in Fig. 5(b). The LoD is subject to scale, application and purpose. LoD might be decided by digitizing regulations (or rules) and accuracy standards published by official mapping agencies in most countries.

Fig. 5.Building models and Level-of-Detail

3.2 Problems of MKP Determination with Surface Fitting

Aerial images and LiDAR data are widely used for building modeling. The entities to determine MKPs are brightness value in images and height in LiDAR data, respectively. LiDAR data provides explicit geometric information while images provide visual and semantic information of the objects. However, fully automatic modeling does not rely on the visual analysis rather requires systematic organizing data to extract geometric information for reconstructing objects. Regardless data source, MKPs should depict distinct, recognizable and representative features of the objects, and must be unique and consistent. However, the results of the extracted MKPs depend on quality of the source data, methodology, implemented algorithms, and required LoD. Several methods are available to extract MKPs such as Douglas-Peuker algorithm, Harris corner detector, Moravec operator, Förstner’s interest operator, chain code and so on (Mikolajczyk and Schmid, 2004; Szeliski, 2011). The important properties of the MKP detection algorithm are scale and rotation invariant. However, none of the algorithms does provide perfect results in most cases.

The main goal of this study is to propose methodology to correct erroneously determined MKPs. Therefore, evaluation of the algorithms extracting MKP is out of scope. In other words, this study focuses on refinement of the MKPs, i.e., relocation of the wrong MKPs to the correct positions. If there are n segmented patches in a building, all possible number of the equation sets to determine a model key points is to be nC3. The key is to find an equation set that provide correct MKPs. In the case of the simple gable roof as shown in Fig. 6, there are six facets (i.e., two roofs and four side walls) that are derived by surface fitting. In this case, numbers of the solutions to determine Model Key Line (MKL) or MKP are 6C2 = 15 or 6C3 = 20, respectively.

Fig. 6.Problematic determination of model key features

However, only one set of the surface equations is acceptable as shown in Fig. 6(a) while the other sets of the equations provide wrong solutions as shown in Fig. 6(b). As a way to solve such problem, Avrahami et al. (2008), Li and Guan (2011), and Sampath and Shan (2008) have suggested utilizing adjacent matrix for building extraction and roof reconstruction. To determine the model key features automatically, adjacent relationships among the facets should be known beforehand. However, to create adjacent matrix is costly process and the reliability might not guaranteed in some cases. This study, on the other hand, adjustment method of MKPs is the main issue regardless extracting MKP methods.

3.3 Determination of MKPs

In this study chain method was applied to extract MKPs as well as to segment clouds. The chain codes were assigned along boundaries of the segmented pathces. However, regularization of the boundaries (i.e., refinement of the boundaries) is required otherwise unnecessary MKPs are extracted. The regularization was easily performed by counting run-length and analyzing pattern of the chain codes. If the run-length is shorter than a certain distance, the points are not taken as MKPs (Fig. 7(a)). However, chain code is not rotation invariant. In order for chain code to be rotation invariant, differential chain codes were computed. Arbitrarily rotated straight line segments have two dominant chain codes and the differences of the codes are the same due to the characteristics of the chain code (i.e., difference of the codes is “4” for opposite pixel with respect to the center). Fig. 7(b) shows the MKPs of the rotated patch. Fig. 7(c) shows detected MPKs of the segmented patch from real data (Engineering building in University of Calgary).

Fig. 7.Detection of model key points

 

4. Correction of Model Key Points

Intersections of the segmented patches after surface fitting for each patch could determine the model key features (i.e., boundaries or corner points of the roofs). In most cases, however, delineation of the segmented patches does not represent the actual roof patches due to imperfect segmentation and interpolation, boundary problem, and other objects around buildings such as vegetation. Incorrectly determined MKPs influence both planimetric and vertical accuracy as shown in Fig. 8.

Fig. 8.Influence of MKPs on object modeling

Surface fitting was performed using segmented data. Planimetric coordinates (X, Y) of the MKPs were substitute into the fitting functions, then the elevation coordinates Z at the MKPs were compared with the original Z coordinates. If the differences were larger than a permissible tolerance (e.g., topographic mapping accuracy standard), the original Z coordinates were replaced by Z coordinates from the fitting functions. However, (X, Y) coordinates of the original MKPs might be incorrectly determined. Therefore, (X, Y) coordinates also should be corrected. In this case, (X, Y) coordinates of the original MKPs were replaced by (X, Y) coordinates that have correct Z coordinates computed from the fitting functions. In order words, locations of the MKPs were determined by mutual refinement of (X, Y) coordinates and Z coordinates.

When fitting is performed with segmented points, it is probably better to exclude some points around edges of the patches because there are points influenced by interpolation and/or trees around the buildings. Therefore, two cases, i.e., using all points in the segmented patches and part of the points in the patches are considered. The points in the areas between roof boundaries and two times of the GSD inside from the boundaries are not used. Results from both cases (i.e., using all points and using central points were compared) were analyzed.

In order to correct erroneously determined MKPs, surface fitting with segmented points was performed. The errors of the MKPs occur (X, Y) coordinates as well.

The original point clouds without gridding were utilized to correct the planimetric coordinates of the MKPs. The correct coordinates were determined by searching the nearest point from the original data. The proposed method could be applicable to both irregular and gridded point cloud data. Fig. 9 shows flowchart with algorithm of the proposed method.

Fig. 9.Flowchart and algorithm to correct MPKs

 

5. Total Least-Squares Solution

Total least-squares solution (TLSS) is a least-squares method of errors-in-variables (EIV) in which observational errors on both dependent and independent variables are taken into account. All (X, Y, Z) coordinates in the point cloud data obtained from ALS system are subject to error. Therefore, object modeling and accuracy evaluation by TLSS is theoretically more reasonable than regular least-squares solution (LESS) because random errors are included in both independent (i.e., X and Y coordinates) and dependent (i.e., Z coordinates) variables. The mathematical model of TLSS is based on Gauss-Helmert Model (GHM) represented as Equation (1) (Schaffrin and Snow, 2010).

where y is observation vector, A is an desig n matrix, eA is an error vector of the independent variables, ξ is unknown parameter vector, and e is error vector of the observations with number of the observations (n) and number of the unknown parameters (m). Matrix A is full rank matrix, i.e., rank(A) = m < n. The expectations and dispersions are E{e} = 0, E{eN} = 0, , and respectively. Assuming there is no correlation between independent and dependent variables, the covariance is to be C{e, eN} = 0.

Hence, Equation (1) could be rewritten as Equation (2).

where ⊗ represents Kronecker product. Therefore, the target function with condition to minimize (eTe+eATeA) is provided by Equation (3).

where λ is Lagrange multiplier. The observation equation of TLSS is represented by Equation (4), and the unknown parameters are computed by Equation (5). Especially, modeling of the curved surface was evaluated by regular LESS and TLSS, and the results were analyzed.

 

6. Results and Analysis

6.1 Test LiDAR Data

Simulated and real airborne LiDAR data were used. The simulation data include polyhedron buildings with gable, pyramid, and dorm (i.e., hemisphere) roofs. The real data is Engineering complex in University of Calgary campus with complicated shapes of the buildings. The average Ground Sampling Distances (GSDs) of the simulated and real data are 0.25m and 0.6m, respectively (Fig. 10). The advantages of using simulation data are to crate various roof types that might not be available in the real data. In addition, quality evaluation can be performed with the simulation data since the ground truths are exactly known.

Fig. 10.Airborne LiDAR data

However, the simulation data may not reflect the real world environment such as vegetation, cars near the buildings, and other objects on the ground. It is obvious that the MKPs of the simulation data are known. The ground coordinates of the MKPs of the real data were available from digital building models. Fig. 11 and Fig. 12 show segmented patches and model key points result from the adaptive 3D chain code method that was introduced in Park (2012).

Fig. 11.Segmented patches

Fig. 12.Model key points

6.2 Correction of Model Key Points and Accuracy

Results of the corrected MKPs and their accuracy are listed in Table 1. The correction method was applied to both using gridded and non-gridded original data for each segmented patch. The coordinates of the MKPs were extracted from gridded data. It is noticed that only Z coordinates were corrected based on the gridded data while both (X, Y) and Z coordinates were corrected based nongridded data. Parts of the results are displayed in the Tables due to the page limit. Fig. 13 shows errors (or residuals) of the Z coordinates for all dataset compared with the true data. The overall RMSEs for all cases, (i.e., gridded and non-gridded, and using all data and using central data) are presented in the bottom line of each Table, and visualized in Fig. 14 for better analysis.

Table 1.Correction results (unit: m)

Fig. 13.Residuals of Z coordinates of MKPs

Fig. 14.Overall RMSE of Z coordinates of MKPs

In most case, to use the points in the central regions of the patches provided higher accuracy as expected. Especially, the overall height accuracies were noticeably improved: RMSEs from 4.043m to 0.114m, from 1.198m to 0.010m, and from 1.333m to 0.417m for pyramid with gable roof and dome with gable roof in simulation data, and real data, respectively. Considering GSD of the dataset and practical accuracy of the ALS system, the RMSEs after correcting erroneous MKPs are impressive. Therefore, proposed method is quite robust. Another noticeable fact is that (X, Y) coordinates of the extracted MKPs are quite accurate while Z coordinates have larger errors as shown in Table 1. This is caused from interpolation and imperfect segmentation.

6.3 Curved Roof Surface Modeling

Planar surfaces are easily modeling with 3D plane functions using Equation (6) while the curved surfaces (e.g., hemisphere, half-cylinder, or other non-planar surfaces) require linearization of the non-linear functions with initial approximations and iteration process. Smoothly curved surfaces do not have distinct MKFs to represent the surfaces. Therefore, modeling of such surfaces could be performed by fitting with an appropriate function such as sphere, halfcylinder, polynomials, or harmonic function that produces minimum RMSE as shown in Fig. 15. Equations (7) and (8) are hemisphere functions for a dome roof modeling. For the automatic curved roof modeling, dome and arch roofs could be identified by analyzing 3D slope vectors. Dome has evenly distributed slope vectors in all directions while arch roof has two dominant directions as shown in Fig. 16 (Lee and Lee, 2010).

Fig. 15.Curved roofs and modeling functions

Fig. 16.3D slope vectors of curved surface

where (X, Y, Z) are coordinates of the point clouds, and a, b, c, R are coefficients or parameters of the surface functions to be determined by least-squares fitting. ao,bo,co, and Ro denote the initial approximations of the hemisphere function. Especially, TLSS applied to the dome roof and compared with results from the regular LESS. Variance component (σo) and Mean-Square-Error (MSE) of the parameters in LESS are computed by Equation (9) and (10), respectively, while σo and MSE of the parameters in TLSS are computed by Equation (11) and (12), respectively.

where v denotes residual. n and m are number of observations and parameters, respectively (Yoo et al., 2011). Table 2 and Fig. 16 show modeling results of a dome roof. The results verify that the TLSS provides reliable parameters compared with the regular LESS based on the various components and MSE of the parameters: Various components of the regular LESS and TLSS are 0.2538m and 0.0255m, respectively.

Table 2.Modeling parameters and accuracy (unit: m)

Fig. 17(a) shows the segmented patches. A curved surface is segmented into several patches because the similar slopes are coded as the same code number. Extracted MKPs from the segmented patches are displayed in Fig. 17(b). Modeling by TLSS fitting with the MKPs without any correction is shown in Fig. 17(c). Finally, Fig. 17(d) shows modeling results using original LiDAR points that are closest to each MKP.

Fig. 17.Dome roof modeling

 

7. Concluding Remarks

ALS technology has advanced in recent years due to the demand from the various applications. Positional accuracy of the LiDAR data fulfills the requirement standards for largescale topographic mapping. However, to develop automatic modeling method is challenging task. Especially, modeling man-made structures is a central issue in LiDAR data processing. Several sophisticated steps are involved in the processing and each step affects the quality of the modeling. The common problems in object modeling are: Segmentation is a crucial step in LiDAR data processing to analyze geometric characteristics of the surface. Automatic extraction of the accurate MKPs from segmented patches is essential for object modeling. Incompleteness of both segmenting patches and extracting MKPs leads inaccurate modeling.

This paper proposed a scheme to correct erroneously determined MKPs to improve accuracy of the building modeling. Following conclusions were drawn throughout this study:

Chain code used as a shape descriptor could be applicable for segmenting and extracting MKPs by extending to 3D space with adaptive and rotation invariant approach. Even though MKPs should be well-defined points to reconstruct objects, to extract accurate MKPs from the segmented patches is challenging task. Therefore, refinement of the erroneous MKPs is required not only for improving modeling accuracy but also quality assessment. Extracted MKPs were evaluated by comparing with true data (i.e., accurate building model data). Since the proposed method are not based on surface fitting, coordinates of the MKPs were compared with also fitting results. In most cases, errors occur along the segmented patch boundaries. Therefore, both fitting cases (i.e., using all data, and using central data excluding data around boundaries) were compared. It is obviously shown that to use the central data provides better accuracy. Almost none of the LiDAR point could represent MKP except the laser pulses hit the exact location of corners of the buildings. However, this is rare case. To evaluate planimetric accuracy of the extracted MKPs, RMSEs were computed using true data. In addition, the closest LiDAR points to the MKPs were selected and the accuracy was evaluated. Accuracy of the extracted MKPs compared with true data:

- Simulation data A: planimetric RMSE = 0.229m, vertical RMSE = 4.043m- Simulation data B: planimetric RMSE = 0.070m, vertical RMSE = 1.198m- Real data: planimetric RMSE = 0.092m, vertical RMSE = 1.333m- The planimetric accuracy is quite acceptable, however, since the vertical accuracy is lower, correction of Z coordinates could be done by replace Z coordinates of the MKPs by fitting results with using central data whose vertical RMSEs are 0.114m, 0.015m, and 0.417m for simulation data A, data B, and real data, respectively. Modeling for different shapes of the curved roof surfaces (i.e., dome and arch) was performed by analyzing 3D slope vectors. However, well-defined MKPs of the curved surfaces could not be detected. Hence MKPs from several patches were used for surface fitting. Both regular LESS and TLSS provided almost the same modeling estimations, however, TLSS might be reasonable approach because all of (X, Y, Z) coordinates of the LiDAR data have errors, i.e., errors in all variables.