Terrestrial and Aerial Laser Scanning Data Integration Using Wavelet Analysis for the Purpose of 3D Building Modeling

1.1. Level of Detail

When creating 3D city models of big cities some methods of 3D building generalization are required [2]. As described in [2], the first and most popular standard, called CityGML, was developed in Germany and categorized in Levels of Details (LoDs) [3]. This standard defines LoDs for buildings models from LoDs 1 to 4. LoD 1 contains coarse building structures, LoD 2 is coarse building blocks with roof models, LoD 3 is filled with detailed building models of building facade, walls, windows, doors, etc., and LoD 4 represents highly detailed architecture building models of LoD 3 but with indoor elements (stairs, installations, detailed elements, figures, sculptures, etc.) [2].

To generate a 3D model of the city with level of details corresponding to LoD 3 and LoD 4, specific conditions have to be fulfilled. The most important are those related to keeping details (generalization) and precision of these models.

One of the main applications of Terrestrial Laser Scanning (TLS) and Aerial Laser Scanning (ALS, sometimes also called LiDAR) is 3D city modeling. Most Digital Surface Models (DSMs) and 3D city models are generated from aerial or satellite images and LiDAR. These models are information rich about the shape of the building, but they usually do not contain many details concerning the building facade. Therefore, such models are textured with additional terrestrial or aerial images.

Currently, a lot of 3D models are being built on the basis of LiDAR data. Aerial Laser Scanning systems can acquire over 40 points per square meter and register multiple echos, which enables accurate terrain and surface reconstruction. However, ALS data is very often incomplete. It means there is a lack of points in some specific areas (building facade, footings, building elements located under tree cover, occlusions, etc., Figure 1A). Therefore, 3D building models could be generated on the basis of terrestrial images [4], but generating models in very narrow streets using terrestrial images can be impossible. In most situations TLS can be the best source of data as it yields very dense and accurate point clouds (1–5 cm) complete with the information that is not acquired by ALS sensors (building interior, facade, etc. Figure 1B). Therefore, complete and accurate 3D modeling of buildings is possible by merging TLS and ALS data.

Additionally there is a question of how to merge two such different types of data: point clouds from TLS and ALS. Nowadays there is still a lack of evaluation systems that can do this task.

The objective of the proposed method is the modification and completion of ALS data set with the use of TLS data. We will show that our method increases the accuracy and density of ALS point clouds.

1.2. Related Work

The issue of TLS and ALS data integration is an active research area [5]. The input information can be divided into two groups: point clouds and imagery from TLS (plus mobile mapping), point clouds and imagery from ALS and 2D maps as well. Then, we can choose the specific data in specific sequences we want to combine. Many approaches include a combination of point clouds, aerial imagery and 2D maps or vector data for building extraction and reconstruction or 3D city modeling [6–8]. Compared with photogrammetric imagery, point clouds of LiDAR can give very accurate XYZ coordinates information. Nevertheless, the quality of boundary lines is quite poor. What is more, in some regions in LiDAR point clouds there might be no data. This situation is caused by self-occlusion [9]. Currently, the most popular combination is combining of point clouds with aerial images. The authors of such papers like [5,10–12] and others, find this integration has many advantages, especially in orthorectification or texturing of 3D models. Another possibility is an approach described by Kremer and Hunter in [13].

According to [14] the basic step which allows us to join photogrammetric data is the orientation to common uniform coordinate system. Rönnholm in [15] suggests a few methods of implementation of one common system for geodata coming from different measurement systems:

-

usage of one hybrid tool which, after proper calibration, allows the gathering of data in one coordinate system simultaneously;

-

preparation of orientation image and laser data separately but in an identity coordinate system;

-

relative orientation, so called registration of different data types.

One of the methods, which allows us to make relative orientation for independently acquired data, is the transformation between systems. Transformation is based on reference points. In this case points are common for both data sets. These points could be selected among natural points in the form of roof corners, edges or other elements which make topographic objects distinctive. Target points also can be made from targets or ground signs [16].

A trend of fusion of different geodata types, which is described for instance in [17–22] has been observed for several years in the Polish and international literature. The most important conclusion from the report is a confirmation of hypothesis about the influence of both a priori accuracy of image and laser data and method of searching for tie points.

In all methods problems related to ALS data can be observed as a result of their “stripped” nature. The only way to adjust and to reduce the influence of those strips to accuracy of transformation is to use affine 3D transformation. There are many examples of such solutions, however, the issue of development of cities' 3D models based only on integrated point cloud from laser scanning is not so common.

In [20], a method of combining data from TLS and ALS using terrestrial and aerial scanner georeferences in a common global coordinate system was described. A terrestrial laser scanner was integrated with a low precision GPS receiver (2 m) and an electronic compass with real precision in urban area equal to 6°, which allowed researchers to get point clouds oriented to a global coordinate system with an accuracy not greater than 2 m in height and 3 m in horizontal coordinate. Moreover, during the integration of data coming from different positions, an Iterative Closest Point (ICP) algorithm was used to increase the accuracy of the general georeference. Afterwards, again with the application of the ICP algorithm, oriented TLS data was joined to DTM and 3D models of buildings were generated from ALS point cloud and cadastral data. Then, the general orientation of TLS scans could be used as a basis for an iterative solution to adjust the two sets of data. Results and analysis of the precision of such solution are described in detail in the mentioned publication.

In the case of integration through data orientation based on the direct transformation of sets of points to common coordinate system (usually global system defined by ALS). Coordinates are recalculated using a specific function with parameters based on reference points. These points are selected among so-called tie points or pseudo-homological points. Then the basic problem is to find such points with high accuracy. Depending on the entry data available, we could distinguish general orientation (a technique used to estimate parameters of transformation between two data sets without knowledge of initial reference points) and detailed orientation (when initial parameters of transformation are already known). General orientation is considered as one of the most popular methods of joining 2D or 3D data sets. The main assumption of such an approach is to define the identity and characteristic elements of the subsets joined. A common area of cover is then required. In such an area some characteristic elements can be distinguished such as points, straight lines [23], curves [24], planes, figures, etc. Another method is an algorithm, described by Chua [25], which defines signature of point. A selected point is described indirectly as a vector calculated by other points located in its neighborhood (defined by particular distance). Then, all selected points are compared with the second set in order to find identity points. In most cases, the accuracy of general registration is not enough to be used in the photogrammetric field. That is why detailed methods are used. However, in order to apply them, initial knowledge about general or precise information about reference points or parameters of transformation has to be available. Common iterative methods are used, which are based on minimizing defined values, like transformation error or defined physical parameters, e.g. the distance between points.

2.1. Methods of Characteristic Points Extraction

In this article a new method of defining coordinates of characteristic points using wavelet analysis is described. Initially, selected ALS and TLS points representing roof edges are the basic component of the algorithm.

2.2. ALS Data—Identification of Edge Points by Indirect Method Using Wavelet Analysis

This kind of integration is based on the transformation between high resolution TLS data and low accuracy ALS data. The extraction of a signal characteristic can be done mathematically and analytically via approximation, averaging, examination of signal's parts or transformation—herein Fourier transformation.

Wavelet analysis represents a signal using hierarchical resolution to receive information about the signal that can be presented in different levels of details. It extracts a general approximation but also enables the selection of some details, disturbances or noise. Theoretically, the variability of relative locations of adjacent points representing the roof edge should be exactly the same in both cases. The main aim of this research was to define whether it is possible to “improve” aerial laser scanning data based on data from terrestrial laser scanning.

Wavelets are sets of functions where each of them is extracted from basic function using shifting and scaling. Generally speaking, wavelet analysis can be described as [26,27]:

Variable s is a scale and τ is shift. Then, reverse transformation could be described as:

Using wavelet analysis, decomposition of signal and its lossless reconstruction is possible [26].

2.3. Algorithm Description

The algorithm specifies the main stages of data processing described below:

3.1. Data Sets

The data used in our research were acquired using TLS and ALS systems (Figure 9). The TLS data sets were acquired using a Leica ScanStation2 in the Military University of Technology area, in the city of Warsaw, Poland, in 2010. This is a pulsed HDS scanner with a range of 300 m and high measurement accuracy (position: 6 mm, distance: 4 mm and high scan resolution below 1 mm). The system is equipped with a terrestrial digital camera (6.4 Mpx). The main issues related to terrestrial laser scanning data acquisition are related to range and incidence angle, that determines indirectly how many points in the upper parts of the building (for example roofs) will be collected. The RMSE of the measured points coordinates was within the range of 3 to 8 cm. To compare we used also a FARO Focus 3D phase laser scanner.

The ALS data sets were acquired in 2009. Scanning was conducted using the ALS Lite Mapper 6800i with a 4–6 points per square meter resolution. ALS measurements data usually consist of not only points information from the top surfaces of the structures (buildings' roofs, the canopy of trees) but also from some branches and lower parts of the trees. This feature will be further named as data incompleteness. Also the merging of objects is very common. The accuracy of measurements was determined using check points (RMSE < 0.30 m in altimetry and 0.10–0.50 m in horizontal error).

The test field consists of buildings with different roof types (flat, steep, arched, etc.) and different structures.

3.2. DWT Signal Processing

Figure 10 presents data divided into vertex and edges with descriptions (k—means edge for one of test building).

Below there is an example of edge kW processing. In the first step, the alignment of the signal was conducted. The results of processing using two proposed method are showed in the Figure 11A,B.

In Figure 11A,B, points in blue represent completed edge data and points in gray present original ALS points. In both cases, the number of edge points is the same. The difference is invisible in the linear and non-linear data character. Below we give an example where the X coordinate of a signal is decomposed for selected edge kW of aerial (kW_A) and terrestrial (kW_T) laser scanning aligned with Method 1. Signals kW_A and kW_T were analyzed in different variants (wavelet family, level of decomposition).

3.3. Reconstruction of Signal and IDWT

As an example of mixed signal reconstruction, we chose another, shorter signal kE3 and found the best decomposition with wavelet “sym3” at level l = 6. For instance, to reconstruct a signal, approximation of original ALS signal from level 6 was used along with details of TLS signal. It could be then understood as “making noise” for the ALS signal based on the TLS signal. The observation of the decomposed signal and its wavelet coefficients allowed us to conclude the following possibilities of reconstruction using an A—approximation of ALS, dT—details of TLS, dA—details of ALS:

(a)

aA6 + dT1 + dT 2 + dT 3 + dT 4 + dT 5 + dT 6

(b)

aA 6 + dT 1 + dT 2 + dA3 + dA 4 + dA 5 + dA 6

(c)

aA 6 + dT 1 + dT 2

Figure 16 presents d1–d6 details for signals of the X coordinate from aerial (XA_kE3) and terrestrial (XT_kE3) point clouds representing the kE3 roof edge. Decomposition was made with sym3 wavelet at level 6.

The difference between both signals as it could be noticeable at lowest frequencies was largest for details d1 and d2 (Figure 16). These differences indicate the different density of “points of roof edge”. Signal reconstruction (XAT_kE3_6sym3_a6d12) was done using approximation a6_XA (coarse ALS signal representation) at level 6, details from ALS: dA6, dA5, dA4, dA3 and details from TLS dT2 and dT1:

• XAT_kE3_6sym3_a6d12 = a6_XA + d6_XA + d5_XA + d4_XA + D3_XA + d2_XT + d1_XT

The signal of X coordinate of kE3 edge was decomposed with the same wavelet “sym3” at level 6. The results are depicted in Figure 17.

Differences could be seen not only for d1 and d2 details but also for higher details. Analysis was made for all coordinates X, Y, Z separately. Finally, the signal of a complete kE3 edge was reconstructed taking into account details d1 and d2 from the terrestrial laser scanning data and presented in Figure 18. The points marked with black colour represent data before decomposition, blue colour represent data after decomposition. More detailed views of edge points are presented in consecutive layers in Figure 19.

Similarly, decomposition and reconstruction were performed for all coordinates for all other edges of test building using both first and second methods of ALS data completion. After edge reconstruction, the vertices of the roof were determined as an edges intersection. These vertices were used as a reference, homologous points in further transformations.

5.1. Quality of the Method

Additional analysis was performed to examine the quality of the proposed method of pseudo-homologous point determination. Table 2 presents the list of edges and segments (column: Edge), lengths (column: Reference Measurement T) calculated on the basis of tachometric measurements and absolute values for edge length differences (column: Difference). These differences were calculated between reference measurement and raw ALS edge data (T-ALS_D), ALS edge data after wavelet reconstruction using first and second method of completing the signal (AP_VA I and II) and edge from TLS data (TLS).

The absolute mean difference value is presented in the last row of Table 2. For edge lengths calculated directly from raw ALS point cloud, the mean value is about 20% bigger than the length from data processed using wavelet analysis. In comparison to particular length values the results are similar. Only two shortest edges of the building cannot be processed using described method. In ALS data the biggest difference does not exceed value of 0.645 m, after wavelet analysis, 0.107 m for the same edge, and 0.427 m generally. Compared to the rest presented methods, where the time of data processing was much shorter, the errors of final transformation were much smaller.

The method was implemented with the MATLAB Wavelet Toolbox. The input data came from point processing according to the method described in Sections 2.1.1–2.3.2. Depending on the chosen method, edge processing is manual or automatic. Our method is semi-automatic, because it needs initially prepared input data. If edges are primarily prepared, the only manual step is determination of wavelet analysis parameters.

The computational cost of the proposed method depends on the data preparation level and quality. For example for unfiltered data there is a need to separate point edges from noise. This can be done manually or automatically. The most time-consuming step is the edge identification and the alignment of the signal. With the implemented algorithm, for one edge it takes about 10–20 min to extract, complete and load data to a MATLAB application. Further steps' time depends on the knowledge and an experience of operator and complexity of data and could take a few minutes. All processing and implementations were made using 2.50 GHz Pentium processor with 4.00 GB RAM and 64-bit WOS.