Registration of Airborne LiDAR Point Clouds by Matching the Linear Plane Features of Building Roof Facets

Abstract: This paper presents a new approach for the registration of airborne LiDAR point clouds by finding and matching corresponding linear plane features. Linear plane features are a type of common feature in an urban area and are convenient for obtaining feature parameters from point clouds. Using such linear feature parameters, the 3D rigid body coordination transformation model is adopted to register the point clouds from different trajectories. The approach is composed of three steps. In the first step, an OpenStreetMap-aided method is applied to select simply-structured roof pairs as the corresponding roof facets for the registration. In the second step, the normal vectors of the selected roof facets are calculated and input into an over-determined observation system to estimate the registration parameters. In the third step, the registration is be carried out by using these parameters. A case dataset with a two trajectory point cloud was selected to verify the proposed method. To evaluate the accuracy of the point cloud after registration, 40 checkpoints were manually selected; the results of the evaluation show that the general accuracy is 0.96 m, which is approximately 1.6 times the point cloud resolution. Furthermore, two overlap zones were selected to measure the surface-difference between the two trajectories. According to the analysis results, the average surface-distance is approximately 0.045–0.129 m.


Introduction
Registration is one of essential preprocesses for the data analysis of point clouds acquired by a laser scanner because the data can be obtained by using different platforms e.g., terrestrial or mobile laser scanning from the ground or airborne laser scanning from airplanes.The point clouds captured by different platforms have good complementarity with the data registration.Different scan stations or trajectories lead to misalignment of a point cloud.The co-referenced point cloud is highly referenced in 3D modeling, surface reconstruction, and information extraction [1][2][3].
The most common approach for point cloud registration in the existing research literature involves artificial targets, such as reflective balls.These targets are regarded as corresponding objects in two different data sources and then apply coordinate transformation models, such as the Bursa model and the 3D rigid body model.In this type of method, the 3D rigid body model is usually assumed, and the scale parameter between two datasets is ignored due to the direct measurement of laser scanning.The method is widely used in terrestrial laser scanning, and many software packages support the automatic extraction of the reflective balls from the point cloud and calculate the coordinates of the ball center.The spatial distributions of the artificial targets are important for model parameter calculation and accuracy evaluation [4].
Remote Sens. 2016, 8, 447 2 of 17 However, artificial targets can hardly be deployed in the measurement campaign of mobile and airborne laser scanning for several reasons.First, applying artificial targets within the measurement field carries a high cost in terms of labor and money because the field is normally large.Second, the density of a point cloud is lower than that of terrestrial laser scanning, especially in airborne laser scanning cases.Therefore, it is difficult to identify target balls from an airborne point cloud.
To co-reference the point cloud without using artificial targets, some automatic registration approaches have been proposed by different researchers.Among them, the ICP (Iterative Closest Point) based algorithm [5] is one of the most popular methods.
The ICP algorithm is sometimes applied to the dataset with different reference systems and aligned in a pre-process.In addition, it is often used to match two datasets with overlapping areas and is usually captured by the same type of platform.The pre-alignment operation is applied to reduce the difference between two datasets and is usually conducted manually.Pre-alignment also provides the rough registration results for precise matching.Then, an iterative refinement technique is used by identifying corresponding points and estimating the best relationship between two datasets.During the relationship estimation, an error matrix is established to evaluate the error between corresponding points.
Based on the basic ICP algorithm, some extensions and additional conditions were developed to increase the robustness and convergence, as well as to improve the computational efficiency and performance.For example, RANSAC (random sample consensus) [6] or LMS (least median of squares estimators) [7] were proposed and applied for data registration.
Another category of approaches for point cloud registration is a geometric-feature based method.Point features, line features, plane features, curvature feature, and terrain should be first extracted and then used to co-reference the different sources of the point cloud.The geometric feature-based algorithms are usually performed semi-automatically because the correspondence among features are normally manually identified.For example, Cheng et al. [8] and Zhong et al. [9] used corner points of buildings and proposed semi-automatic method for point cloud registration.[10] proposed a semi-automatic method using feature line and feature polygon to matching the point cloud captured from terrestrial and airborne laser scanning systems under the assumption that the rotation angles between X and Y are zero, which is difficult to satisfy in a real situation.However, this method can be used to register a point cloud with low overlap zones.Li and Guskov [11] proposed a method to detect a set of salient feature points and use a significance vector as an approximate invariant.Then, the feature sets were used to obtain the alignment between two datasets.Similar to [11], Alba et al. [12] and Weinmann et al. [13] used line and plane features extracted from a point cloud to obtain accurate alignment results.Variation in the curvature is a special feature used for point cloud registration [14].The normal vector and the curvature from an unorganized point cloud should be calculated in advance [15].Hansen et al. [16] used the 3D line features derived from both terrestrial and airborne point clouds registration.Then, a two-step registration method is performed.According to their research, this method will be easily affected by noise.A recent publication [17] proposed a registration method based on a hierarchical and coarse-fine strategy to co-reference a point cloud from airborne and vehicle platforms.According to their analysis, this method can achieve an average accuracy of approximately 0.4 m.
Curvature features [18] are other types of features used for registration.The curvature parameter for each point is computed from different datasets, and the initial feature points are selected by maximum local curvature variations.Last, the Hausdorff distance is used to obtain accurate matching points and to estimate the rotation and translation matrix.Some other types of data sources are also used to aid point cloud registration, such as imagery and topographic maps.Dold and Brenner [19] proposed an autonomous method based on planar patches automatically extracted from point clouds and imagery.Wu et al. [20] proposed a method using a topographic map as reference data and transferred a point cloud captured from both terrestrial and airborne laser scanning systems to the local coordinate system.During this research, the building feature lines and control points of topographic map were selected for feature based registration.Supplementary data sources usually imports a local or special reference coordinate system; therefore, this type of method is helpful for further integration and application.[21] proposed a 3D building histrogram to represent the complexity of buildings and uses SD (block distance between the two histograms) and ZLB (z-value of the largest class in the 3D histogram) as the indices to describe each building roof and co-register the multi-layer and single-layer buildings.Yu et al. [22] extracted the semantic features (such as facades, poles, cars) from city environment and used the ICP algorithm to align the point cloud from Google Street View cars.
This paper proposes a novel approach for the registration of airborne point cloud data sets by using linear polygon features of building roofs.In the first step, the methods developed in previous research [23] will be applied to extract buildings and identify buildings with gabbled and hipped roofs in the two datasets, respectively.The buildings identified with the same building footprint on OpenStreetMap (OSM) are, therefore, the corresponding buildings in the subsequent process.Next, the normal vectors are calculated for every simple roof facet of the corresponding buildings.These normal vectors will be used to calculate the transformation model between the two datasets; therefore, the least squares (LS) method will be deployed to solve the observation system with redundant feature-pairs.The rest of the paper is organized as follows.Section 2 introduces the proposed method, including the principles of the registration method in this paper, the flowchart of the proposed method, and the mathematical models.In Section 3, the study area, as well as the experimental results and the accuracy evaluation results, is presented.Section 4 summarizes the proposed method and discusses the limitations of the proposed method.

Point Cloud Registration by Matching Plane Features of Building Roof Facets
To describe the two dataset used in this paper, we discuss two point cloud datasets to be registered as ´X Y Z

¯C1
and ´X Y Z

¯C2
, where the subscript C1 and C2, respectively, are used to separate the datasets.The purpose of this paper is to find a proper transformation parameter to reduce the difference between two datasets.In this paper, the rigid body model is chosen for point cloud transformation, whereas the steps of the other models (e.g., Bursa model) for registration are similar to the rigid body model.The rigid body model uses a 3D shift vector and a 3D rotation matrix to describe the relationship between two datasets.The mathematical model is described as: where p∆x, ∆y, ∆zq is the 3D translation vector, and R pα, β, γq is the 3D rotation matrix that is defined by the axis angles pα, β, γq between C1 and C2.Therefore, computing the proper 3D translation vector and the 3D rotation matrix is the essence of registration, and it determines the accuracy of registration.

Principle and Flowchart of Proposed Method
For a typical coordinate transformation model, 3D points are usually used to calculate the 3D translation vector and the 3D rotation matrix of the rigid model.After selecting several 3D point pairs from different datasets, a least square adjustment processing is applied to compute the parameters.However, laser scanners use a huge number of points to describe a certain object.The accurate point-to-point relationship is difficult to follow.Therefore, in this paper, the linear polygonal features of buildings, instead of 3D point-pairs, are selected to calculate the coordinate transformation parameters.Figure 1 shows the point cloud of a building roof and its corresponding normal vectors.An airborne laser scanner obtains the reflected laser beams to compute the coordinate of surface of terrain objects.In an urban area, the surface points of building are the main components of an urban point cloud.Meanwhile, linear polygonal roof features are the main features of a building; such features can be easily found and the feature parameter can be easily obtained to perform airborne point cloud registration.
Figure 2 illustrates the process of the registration.After the simple buildings (with simply structured roofs) are extracted, the vectors of the roof facets are calculated.The vectors of the corresponding roof facets can then be used to establish an observation equation system to estimate the transformation between the two datasets.The flowchart can be divided into two main parts.An airborne laser scanner obtains the reflected laser beams to compute the coordinate of surface of terrain objects.In an urban area, the surface points of building are the main components of an urban point cloud.Meanwhile, linear polygonal roof features are the main features of a building; such features can be easily found and the feature parameter can be easily obtained to perform airborne point cloud registration.
Figure 2 illustrates the process of the registration.After the simple buildings (with simply structured roofs) are extracted, the vectors of the roof facets are calculated.The vectors of the corresponding roof facets can then be used to establish an observation equation system to estimate the transformation between the two datasets.The flowchart can be divided into two main parts.An airborne laser scanner obtains the reflected laser beams to compute the coordinate of surface of terrain objects.In an urban area, the surface points of building are the main components of an urban point cloud.Meanwhile, linear polygonal roof features are the main features of a building; such features can be easily found and the feature parameter can be easily obtained to perform airborne point cloud registration.
Figure 2 illustrates the process of the registration.After the simple buildings (with simply structured roofs) are extracted, the vectors of the roof facets are calculated.The vectors of the corresponding roof facets can then be used to establish an observation equation system to estimate the transformation between the two datasets.The flowchart can be divided into two main parts.First, the point clouds of the building roofs are extracted from the two LiDAR datasets.This extraction is achieved by using an OSM aided approach presented in [23].In the next step, the point clouds segmentation will be conducted by adapting the method of ridge-based hierarchical decomposition.Specifically, a roof will be selected for the subsequent process when only less than one ridge can be detected.In other words, only simple structured roofs, such as gabble, hipped, or flat roofs, will be selected for the registration of the two datasets.Then, the point clouds of simple structured buildings are automatically selected based on the roofs.
The second part is the registration parameter calculation and will be divided into two steps: the rotation matrix calculation and the 3D translation vector calculation.The detailed steps and functions for rotation matrix and 3D translation vector calculation will be introduced in Section 2.3.

Calculating the Normal Vectors of Building Roof Facets
A plane-fit program is performed after the building roof points are separated.The detailed steps of the plane parameter extraction are shown in Figure 3.
Remote Sens. 2016, 8, 447 5 of 16 First, the point clouds of the building roofs are extracted from the two LiDAR datasets.This extraction is achieved by using an OSM aided approach presented in [23].In the next step, the point clouds segmentation will be conducted by adapting the method of ridge-based hierarchical decomposition.Specifically, a roof will be selected for the subsequent process when only less than one ridge can be detected.In other words, only simple structured roofs, such as gabble, hipped, or flat roofs, will be selected for the registration of the two datasets.Then, the point clouds of simple structured buildings are automatically selected based on the roofs.
The second part is the registration parameter calculation and will be divided into two steps: the rotation matrix calculation and the 3D translation vector calculation.The detailed steps and functions for rotation matrix and 3D translation vector calculation will be introduced in Section 2.3.

Calculating the Normal Vectors of Building Roof Facets
A plane-fit program is performed after the building roof points are separated.The detailed steps of the plane parameter extraction are shown in Figure 3.After the roof points are extracted from the original point cloud, a plane regression algorithm is applied to the roof points and the plane parameters (see A, B, C, and D in Equation ( 2)) are obtained.Several mathematical forms can be used to represent the plane equation.
where (A, B, C) represent the unit normal vector of the plane, and D denotes the distance between origin to the regressed plane.Some of the roof points will satisfy the plane parameters; however, some points will not.To select the unsatisfied points, a threshold is used to evaluate the residual error for each point.If the residual error of a point is larger than the threshold, then this point will be removed and will not be used for the plane parameter calculation in the next regression.After several iterations, if the residual error of the entire point cloud is lower than the threshold, the regression will be stopped, and then the plane parameters (A, B, C and D) will be used for further processing.
For clarity, we use ( ) and ( ) to represent the plane parameters extracted from C1 and C2 dataset.

Model Parameters Estimation Using Common Building Features
The traditional coordinate transformation model usually adopted 3D point-pairs to establish the relationship between the coordinates of the point pair and the transformation model.Then, after After the roof points are extracted from the original point cloud, a plane regression algorithm is applied to the roof points and the plane parameters (see A, B, C, and D in Equation ( 2)) are obtained.Several mathematical forms can be used to represent the plane equation.
where (A, B, C) represent the unit normal vector of the plane, and D denotes the distance between origin to the regressed plane.Some of the roof points will satisfy the plane parameters; however, some points will not.To select the unsatisfied points, a threshold is used to evaluate the residual error for each point.If the residual error of a point is larger than the threshold, then this point will be removed and will not be used for the plane parameter calculation in the next regression.After several iterations, if the residual error of the entire point cloud is lower than the threshold, the regression will be stopped, and then the plane parameters (A, B, C and D) will be used for further processing.
For clarity, we use ´A B C D

¯C1
and ´A B C D

¯C2
to represent the plane parameters extracted from C1 and C2 dataset.

Model Parameters Estimation Using Common Building Features
The traditional coordinate transformation model usually adopted 3D point-pairs to establish the relationship between the coordinates of the point pair and the transformation model.Then, after introducing the initial value of model parameters, a least square adjustment method is performed to calculate the correction value of model parameters.Last, the optimal value of model parameters is obtained.Similar to the traditional method, the parameter calculation method also builds a relationship between the observational data and the model parameters.Then, the least squares method is applied to calculate the model parameters.

Rotation Matrix Calculation
For convenience, suppose P and Q are two ordinary real points in a building; the corresponding points in the first trajectory (C1) are P C1 and Q C1 , whereas in the second trajectory (C2), the corresponding points are P C2 and Q C2 .According to the coordinate transformation model, with some proper model parameters, the coordinates of two point-pairs satisfy the following relationship: Subtracting the above two equations, we obtain: where `∆X PQ , ∆Y PQ , ∆Z PQ ˘are the three components of vector between point P and Q. Equation ( 4) describes the relationship between homonymy vectors in different datasets.Considering that the normal vector of a building plane is also a type of homonymy vector in two datasets, we can replace the vector between two points by the normal vector of the building plane.Thus, we obtain: where ´A B C ¯T C1 is the normal vector of a plane in C1 dataset, and ´A B C ¯T C2 is the normal vector of the corresponding plane in C2 dataset.The normal vectors in C1 and C2 can be calculated by the approach introduced in Section 2.2.Thus, the purpose of this step is to calculate the unknown matrix R pα, β, γq.
According to Equation (5), as soon as three homonymy normal vectors are selected between two datasets, the night elements of rotation matrix R pα, β, γq will be calculated uniquely.Although more than three homonymy normal vectors exist, the rotation matrix can be calculated by the least square (LS) method.Equation (6) gives the solution when abundant homonymy normal vectors were used for calculation: R " invpV 2 where is the normal vector of a building plane in the C2 dataset, whereas is the corresponding normal vector of a building plane in C1. n denotes the number of normal vectors used during the registration.This method avoids the difficulty in calculating the three independent unknown parameters (α, β, γ) and directly computed the elements of matrix R. As R denotes the rotation matrix between reference and objective coordinate system, an additional requirement of det(R) = 1 should be satisfied when computing R by iteration.

3D Translation Calculation
The method introduced in Section 2.3.1 can be used to obtain the rotation relationship between C1 and C2.However, the shift relationship between two datasets has not been determined.Thus, to calculate the 3D shift parameters, suppose C3 is the dataset that was rotated from C1 by function Equation (7).The corresponding plane parameters are Considering that the fourth parameter (D) is the distance from origin to plane, it will not be changed after the rotation.Therefore, before and after rotation, parameter D will remain the same.Then, we obtain: Meanwhile, after rotation, the normal vector of a plane in C1 will be changed in parallel with the corresponding plane in C2.Therefore: According to the transformation model adopted in this paper, after rotation, only the 3D shift parameters are unknown between C3 and C2.Suppose the corresponding point of in the C3 dataset after rotation is P C3 " ´XP Y P Z P ¯C3 .According to the plane parameters, ´A B C D ¯T C3 satisfy: After introducing 3D translation parameter to function Equation (10), we obtain: The following can be obtained: Substituting Equations ( 8) and (10) to Equation (11), we obtain: According to the above equation, as soon as three homonymy normal vectors are selected between the two datasets, the three components of the vector shift can be calculated uniquely.Similarly, with rotation parameter calculation, although more than three homonymy normal vectors are valid, the shift parameters can be calculated by the least square method.¨∆x ∆y ∆z

Experiments
Two cases with two adjacent trajectories of point clouds were used for the registration experiments.The results are demonstrated in this section.

Details of Datasets and the Selected Corresponding Roofs
Two cases were included in this part.Each case contains two trajectories' point cloud.The first case are located in Toronto and were captured in the year 2011.The two trajectories contained 2,727,981 points and 1,928,393 points.The typical objects in the case dataset were buildings, roads, grasses, and trees.The average point distance of the datasets was approximately 0.6 m.A part of the direct overlapping of these two datasets can be seen in Figure 4.There were obvious offsets and a slight rotation between the two datasets.The translation between the two datasets can be estimated by calculating the distances of a number of manually selected corresponding points.The result was approximately 4.1 m on average.
The second case are located in Zhoushan, a city of Eastern China and were captured in the year 2008.Two trajectories contains 356,489 and 784,271 points.Typical objects in this case were buildings, roads, trees, and water.The point density was 2.4 points/m 2 .The case area is near the airport of Zhoushan; therefore, buildings here were very simple.Flat roofs and gable roofs are two main roof types included in the dataset.The difference between two trajectories' point is relative large.A preliminary evaluation of difference of two trajectories point cloud is about 25 m on average.
In this study, both two trajectories' datasets of two area were used to compute the registration parameters.The selection of building roofs and the alignment parameters were listed in Sections 3.2 and 3.3.A comparison between proposed method and the other methods, like ICP and LMS were also listed in Section 3.3.Then, some further analysis, including the accuracy evaluation and error sensitivity evaluation was furnished by the first dataset because of the complexity of this dataset.The results will be given in Sections 3.4 and 3.5.

Selection of Building Roofs for Registration
First of all, buildings have to extracted from the point clouds.In this work, an OSM-aided approach [23] is applied for this task.Considering the fact that the translation of two datasets is normally not that much, the buildings in the point clouds corresponded with the same building footprint will be regarded as corresponding buildings.In the next step, building roofs are segmented using the method presented in the same paper [23].Only those buildings are selected if their roofs contain none or one ridge.In other words, the corresponding buildings are selected for further process, only when their roofs are simply constructed.Then feature regression method were applied to the point cloud of the selected building roofs.

Details of Datasets and the Selected Corresponding Roofs
Two cases were included in this part.Each case contains two trajectories' point cloud.The first case are located in Toronto and were captured in the year 2011.The two trajectories contained 2,727,981 points and 1,928,393 points.The typical objects in the case dataset were buildings, roads, grasses, and trees.The average point distance of the datasets was approximately 0.6 m.A part of the direct overlapping of these two datasets can be seen in Figure 4.There were obvious offsets and a slight rotation between the two datasets.The translation between the two datasets can be estimated by calculating the distances of a number of manually selected corresponding points.The result was approximately 4.1 m on average.The second case are located in Zhoushan, a city of Eastern China and were captured in the year 2008.Two trajectories contains 356,489 and 784,271 points.Typical objects in this case were buildings, roads, trees, and water.The point density was 2.4 points/m 2 .The case area is near the airport of Zhoushan; therefore, buildings here were very simple.Flat roofs and gable roofs are two main roof types included in the dataset.The difference between two trajectories' point is relative large.A preliminary evaluation of difference of two trajectories point cloud is about 25 m on average.In total, six buildings roofs of the first case were selected for the registration because their roofs are simple structures.In Figure 4, the polygons denote the locations of the selected buildings and the color denotes the elevation of point cloud.
According to previous studies, there are over 33 types of roofs exists in the world, such as flat, gable, hip, L-joint, T-joint, etc. [24].In this paper, three simple types of roofs were adopted: flat, slope, and gabled roof (see Figure 5).These three roof types are very common and typical in urban areas and is easy to fit the plane parameters from point cloud.Moreover, to avoid an ill-conditioned matrix in the observation system when using the least square method, different types of roofs with more points are recommended for use.
In total, six buildings roofs of the first case were selected for the registration because their roofs are simple structures.In Figure 4, the polygons denote the locations of the selected buildings and the color denotes the elevation of point cloud.
According to previous studies, there are over 33 types of roofs exists in the world, such as flat, gable, hip, L-joint, T-joint, etc. [24].In this paper, three simple types of roofs were adopted: flat, slope, and gabled roof (see Figure 5).These three roof types are very common and typical in urban areas and is easy to fit the plane parameters from point cloud.Moreover, to avoid an ill-conditioned matrix in the observation system when using the least square method, different types of roofs with more points are recommended for use.In the next step, the normal vectors of these roof facets were calculated and used to estimate the transformation between the two datasets.The results are shown in 3.3.

Estimation of Transformation Parameters
In this subsection, several buildings of each dataset were selected from the original point cloud and the normal vectors were computed.Then the transformation parameters were then calculated by normal vectors.As a comparison, two methods, ICP and LS, were selected to verify the validity of parameters.The source code of ICP were downloaded from https://github.com/charlienash/nricp.During the LS processing, nine homonymous points for each dataset were manually selected to compute the parameters.In the next step, the normal vectors of these roof facets were calculated and used to estimate the transformation between the two datasets.The results are shown in 3.3.

Estimation of Transformation Parameters
In this subsection, several buildings of each dataset were selected from the original point cloud and the normal vectors were computed.Then the transformation parameters were then calculated by normal vectors.As a comparison, two methods, ICP and LS, were selected to verify the validity of parameters.The source code of ICP were downloaded from https://github.com/charlienash/nricp.During the LS processing, nine homonymous points for each dataset were manually selected to compute the parameters.

Parameters comparison for the first case
Six buildings with 10 different roofs were selected for the first case area.The transformation parameters were computed and listed in Table 1.

Parameters comparison for the second case
Thirteen buildings with 15 different roofs were selected for the second case.Normal vector of planes were calculated and introduced into the proposed method; the transformation parameters, both proposed method, ICP method, and LS method were computed and listed as following.As the difference between two trajectories is obvious, a preprocessing to roughly co-register the two trajectories' data is performed for ICP and LS methods.The results were listed in Table 2. Compared with the above parameters, the rotation matrix of proposed method is quite similar with the ICP and LS method.However the shift parameter is quite different because the difference between two trajectories is quite large.

Accuracy Evaluation
In order to evaluate the accuracy for proposed method, the first case was selected because the buildings in first case is relative more complex.Here, checkpoints as well as the two overlap buildings were selected to verify the difference after transformation.The horizontal distance of checkpoints were computed to evaluate the horizontal accuracy, while the overlap difference were computed to evaluate the vertical accuracy.

Accuracy Evaluation by Checkpoints
To evaluate the horizontal accuracy of registration, 40 checkpoint-pairs were selected, and the checkpoint distances were measured.The checkpoint-pairs were manually chosen in 2D point cloud view and most of them were building corners.The horizontal distance of each checkpoint pair was used as the index to evaluate the accuracy after the registration.Some statistical values, such as minimum, maximum, and average, were also computed.Table 1 shows that the horizontal distance ranged between 0.195-1.628m, with an average distance of 0.964 m, which was approximately 1.6 times the average resolution of the raw point cloud (0.60 m).As a comparison, the same values from point cloud which were transferred by transmission parameters calculated by ICP method and LS method were also computed and shown in Table 3. Results show that the average, minimum and maximum value between checkpoint-pair by proposed method is greater than ICP method and LS method.Different trajectories' point clouds are usually captured from different points of view and overlap zones exist from adjacent trajectories' point clouds.Here two overlapping buildings were selected to evaluate the residuals after registration.The purpose of overlap zone evaluation is to find the registration results in the overlap area by measuring the overlay difference, which also means the The graphical result for the second overlap area is shown in Figure 7.The maximum distance of the case area was approximately 1.100 m, and the average distance was 0.129 m.The graphical result for the second overlap area is shown in Figure 7.The maximum distance of the case area was approximately 1.100 m, and the average distance was 0.129 m.Then, eleven levels were used to divide the points to analyze the distribution of different residuals.The number of points, percentage and cumulative percentage were then calculated; the results are listed in Table 5.According to the table, over 90% of the point-surface distances were lower than 0.30 m, demonstrating good quality in the results.Then, eleven levels were used to divide the points to analyze the distribution of different residuals.The number of points, percentage and cumulative percentage were then calculated; the results are listed in Table 5.According to the table, over 90% of the point-surface distances were lower than 0.30 m, demonstrating good quality in the results.

Error Sensitivity Evaluation of Proposed Method
Section 3.3 demonstrate that the proposed method can be used to compute the transformation parameters.However, the roof features regression processing will always be affected by the sensor error and noise.Therefore, this subsection is designed to evaluate the effect of random errors and demonstrate the stability for the proposed method.
In order to compare the transformation parameters calculated by proposed method with real parameters.We temporarily computed the C3 data which was transformed by C1 with the parameter, p∆x, ∆y, ∆zq true and R true : p∆x, ∆y, ∆zq true " 3748.245, 1569.256, 12.235 The pα, β, γq were set with the value of 1.2 ˝, 2.2 ˝, 3.2 ˝degrees separately.Then, according to the calculation equation, the true value of the rotation matrix is: The p∆x, ∆y, ∆zq true and R true will be treated as the true transformation parameters between C1 and C3.
The transformation parameters were, firstly, computed according to the proposed method, and then the difference between computed parameters and true value is calculated (See first line of Table 6).Then, four levels of random errors were added to each point of C3 to testify the effect of random errors for proposed method.Then, transformation parameters, as well as the difference between computed parameters and real values were calculated again (See line 2-4 of Table 6).The last two columns of Table 6 are calculated by the true value minus the computed transformation parameters in this experiment.We can find that once random errors were added to the original point cloud, the sensitivity of the rotation matrix is obvious because the components' differences of rotation matrix decreased from about 1 ˆ10 ´12 to 1 ˆ10 ´6 rapidly.However, as the random error increased from 0.001 m to 0.100 m, the components difference decreased slowly (from about 1 ˆ10 ´6 to 1 ˆ10 ´4).From the last column of above table, the shift parameter difference is relevant with the scope of the random errors.Once the scope of random error is relatively small, the difference of the 3D shift parameter is low.

Visualization of the Matched Point Cloud
Applying the transformation parameters distributed in Section 3.2, the point cloud of trajectory 1 was transferred.Figure 8 shows the results of a complex building in the case area, which demonstrates that the proposed method achieved good results.The two colors of the figure denote the trajectory number: the orange color denotes the first trajectory, and the green color denotes the second trajectory.
Applying the transformation parameters distributed in Section 3.2, the point cloud of trajectory 1 was transferred.Figure 8 shows the results of a complex building in the case area, which demonstrates that the proposed method achieved good results.The two colors of the figure denote the trajectory number: the orange color denotes the first trajectory, and the green color denotes the second trajectory.

Discussion
Although the mathematical calculations of the proposed method are simple, two strict requirements are required when applying this approach.The first requirement is that the direction of normal vectors of plane-pairs should be the same in both datasets.The normal plane regression method will provide a normal vector ( ) and the constant value (D) as the plane parameters.
These parameters will satisfy the plane equation AX + BY + CZ = D; however, the negative form of the normal vector (− − − ) and the negative constant value ( − ) also satisfy the plane equation.We must choose the proper form of the plane parameters and establish suitable plane-pairs.In this paper, an additional operation was applied to the regressed plane parameters.An angle between the normal vector and the Z axis was used to evaluate the status of normal vector.If the angle was over 90 degrees, then we used the negative form of the normal vector in the relevant equations.
Another requirement to use the proposed method is the selection of a variety of normal vectors, i.e., the building planes used to calculate the model parameters should point in different directions.More than three different directions are recommended.In this paper, we selected three types of building planes, and among six normal vectors, three different directions were used.This will prevent the matrix ( ) from becoming an ill-conditioned matrix and lead to a successful solution.
Meanwhile, as a rotation matrix, an additional requirement, det(R) = 1, should be satisfied during the calculation.

Discussion
Although the mathematical calculations of the proposed method are simple, two strict requirements are required when applying this approach.The first requirement is that the direction of normal vectors of plane-pairs should be the same in both datasets.The normal plane regression method will provide a normal vector ´A B C ¯and the constant value (D) as the plane parameters.These parameters will satisfy the plane equation AX + BY + CZ = D; however, the negative form of the normal vector ´´A ´B ´C ¯and the negative constant value (´D) also satisfy the plane equation.We must choose the proper form of the plane parameters and establish suitable plane-pairs.In this paper, an additional operation was applied to the regressed plane parameters.An angle between the normal vector and the Z axis was used to evaluate the status of normal vector.If the angle was over 90 degrees, then we used the negative form of the normal vector in the relevant equations.
Another requirement to use the proposed method is the selection of a variety of normal vectors, i.e., the building planes used to calculate the model parameters should point in different directions.More than three different directions are recommended.In this paper, we selected three types of building planes, and among six normal vectors, three different directions were used.This will prevent the matrix pV 2 1 V 2 q from becoming an ill-conditioned matrix and lead to a successful solution.Meanwhile, as a rotation matrix, an additional requirement, det(R) = 1, should be satisfied during the calculation.
For most registration research, control points are the common features used to compute the registration parameters.In this paper, plane roof features were used instead of control points, and good results were achieved.However, the detailed principle of adjustment by roof features still requires a mathematical discussion.Moreover, as the real transformation parameters are unknown, the estimated parameters by control points or control planes are approaching the real values.Moreover, the distribution of the selected roof features is important for the calculation of the model parameters.The selection of more features or the altering of other features will lead to different results.Therefore, we suggest that the user selects roof features by equal distribution.

Conclusions
This paper proposed a method using the linear plane of buildings as the sources of feature-based data registration.A three-step calculation method of the model parameters was introduced.The first step of the method is to select the simple-structured buildings with the assistance of OpenStreetMap (OSM).Then, the normal vectors of the selected building roofs are computed using a fitting program and then are introduced into an observation system to obtain the proper translation and rotation parameters.Last, the rotation and shift parameters of the 3D rigid body model are applied to the point cloud to be matched.Two trajectories' point clouds were used to verify the validity and accuracy of proposed method.According to the experimental results, the proposed method achieved good results and can be used to co-reference the adjacent point cloud.
The novelty of this paper is that the linear features of the corresponding roof facets were used for the estimation of the registration parameters instead of the corresponding points.This process avoids the large ambiguity when selecting corresponding points manually and reduces the cost when using artificial objects as the corresponding points.

Figure 1 .
Figure 1.The normal vector of plane roof facets on a building.

Figure 2 .
Figure 2. The flowchart of the proposed method.

Figure 1 .
Figure 1.The normal vector of plane roof facets on a building.

Figure 1 .
Figure 1.The normal vector of plane roof facets on a building.

Figure 2 .
Figure 2. The flowchart of the proposed method.Figure 2. The flowchart of the proposed method.

Figure 2 .
Figure 2. The flowchart of the proposed method.Figure 2. The flowchart of the proposed method.

Figure 4 .
Figure 4. Spatial distributions of the selected buildings.(a) Selected buildings of first case; and (b) selected buildings for the second case.

Figure 4 .
Figure 4. Spatial distributions of the selected buildings.(a) Selected buildings of first case; and (b) selected buildings for the second case.

Figure 8 .
Figure 8. Visualization of the local building's point cloud after registration.

Figure 8 .
Figure 8. Visualization of the local building's point cloud after registration.

Table 1 .
Comparison of transformation parameters for the first case.

Table 2 .
Comparison of transformation parameters for the second case.

Table 3 .
Accuracy of the checkpoints (m).

Table 5 .
Number of points in different layers of overlap zone-2.

Table 6 .
Rotation and shift parameters difference after random error added.