Abstract
Nowadays, terrestrial Lidar scans can cover rather a large area; the point densities are strongly varied because of the line-of-sight measurement principle in potential overlaps with scans taken from different viewpoints. Most of the traditional methods focus on registration algorithm and ignore searching model. Sometimes the traditional methods are directly used to align two point clouds; a large critically unsolved problem of the large biases will be created in areas distant from the overlaps while the local overlaps are often aligned well. So a novel measurement axis searching model (MASM) has been proposed in this paper. The method includes four steps: the principal axis fitting, the measurement axis generation, low-high-precision search, and result generation. The principal axis gives an orientation to the point cloud; the search scope is limited by the measurement axis. The point cloud orientation can be adjusted gradually until the achievement of the global optimum using low- and high-precision search. We perform some experiments with simulated point clouds and real terrestrial laser scans. The results of simulated point clouds have shown the processing steps of our method, and the results of real terrestrial laser scans have shown the sensitivity of the approach with respect to the indoor and outdoor scenes.
1. Introduction
Laser scanners are standard 3D data acquisition devices in a wide range of applications like airborne, mobile, static terrestrial laser scanners, and so forth. In this paper, we deal with LiDAR point clouds of static terrestrial laser scanners (TLS). Our research mainly focuses on search matching for two point clouds registration.
The point cloud registration methods are divided into three groups. First group is named as the ICP (Iterative Closest Point) algorithm [1, 2] and its variants, namely, ICP variants (Yang and Medioni [2] and Grant et al. [3] and Barnea and Filin [4]). This method requires that the two point-could densities in the overlap region must be same. The second group belongs to the approach based on probability and statistical modeling. Montesano et al. [5] presented a registration algorithm based on the same proposition called the probabilistic iterative correspondence method (pIC). Magnusson et al. [6–9] described the NDT (Normal Distribution Transform) method. Granger and Pennec [10] formulated registration process as a general maximum-likelihood (ML) estimation of the transformation and the matches. Segal et al. [11] proposed Generalized-ICP, which unifies the ICP formulations for various error metrics such as point-to-point, point-to-plane, and plane-to-plane. The last group belongs to the methods based on feature vector with multidimensional attributes. The feature vector describes a point and its associated local surface geometry such as spin images [12] and fast point-feature histograms (FPFH) [13].
Nowadays, the terrestrial laser scanner can cover rather a large area and the point densities vary in potential overlaps, although sometimes only a few valid points appear in potential overlaps. The above-mentioned algorithms focus mainly on registration features and ignore the point cloud shape. If such algorithms are directly employed to align laser scans, a large bias will usually be created in regions distant from the overlap area; however, the alignment of local overlaps will be maintained well at the same time. Nüchter et al. [14] computed octree representations from laser scans and aligned them to get an initial estimation for the subsequent ICP scan matching. Franaszek et al. [15] describe a laser scan registration method using sphere targets. However, the placement of the targets is not convenient in some applications.
The failure of terrestrial laser scans registration does not owe to the registration algorithm itself, rather it arises due to falling into a local minimum. Therefore, we propose a search measurement axis searching model (MASM) to register two laser scans. The registration result is corrected and updated continuously through principal axis fitting, measurement axis generation, low-precision search, high-precision search, and optimal result generation. The goal of MASM is to search a better registration position between the two point clouds and avoid the local minimum to keep global geometrical consistency of the registered point clouds. Figure 1 shows the workflow of the proposed algorithm.
This paper is organized as follows: Section 1 describes the basic principle and process of the algorithm proposed. Section 2 examines the sensitivity of the approach with respect to indoor and outdoor scenes. Conclusions and future research directions are discussed in Section 3.
2. Registration Method
2.1. Principal Axis Fitting
Two simulated point clouds (x and y) are used here to describe the overall process of the proposed algorithm. Suppose that the total numbers of points in both point clouds are same covering the same area range. Each point cloud includes 4 points, designated as x-1 to x-4 and y-1 to y-4 forming a total of four pairs of points. These points are illustrated diagrammatically in Figure 2. Here, x is called the dynamic point cloud, and y is called the static point cloud. In this algorithm, all the rotational calculations for point clouds are performed using the quaternion’s method [16, 17].
Firstly, the three concepts involved are described as follows.
Geometric Center. The average value of the coordinates of a point cloud defines the geometric center.
Principal Axis. The principal axis of a point cloud is the line connecting the geometric center to the most distant point. The principal axis represents the orientation of the point cloud.
Geometric shape of a point cloud is described using geometric center and principal axis.
The principal axis shows the initial orientation of the point cloud. When the principal axis fitting process starts, the geometric centers of both point clouds must coincide with each other as depicted in Figure 3. Next, the point cloud x is rotated about the geometric center until the principal axes of x and y coincide with each other, as shown in Figure 4. The result of the principal axis fitting process can act as the initial estimation of the measurement axis generation which needs to be calculated in the later stage.
In general, the principal axis fitting can provide a good initial estimation, but sometimes not. It is not a big deal, just adding more iteration needed in the subsequent measurement axis generation step.
2.2. Measurement Axis Generation
In theory, x has an infinite number of positions with respect to y, which signifies that the search matching is rather time-consuming. The different measurement axes which are produced in each iteration tend to limit the search scope.
In this work, five measurement axes of y, numbered from 0 to 4, are created by the first iterative calculation as shown in Figure 5. Axis 0 is referred to as the “center” axis, and it coincides with the principal axis of y. The other four measurement axes are evenly distributed around axis 0 and are referred to as “side” axes. α is termed as the shrink angle between the side axis and the center axis, having no fixed value, but it will keep changing according to the actual situation. In the first iteration, initially the principal axis of x is fitted with axis 0 of y while x is rotated about axis 0 of y, and the error is computed between x and y. After x has traversed all the five measurement axes, five errors are obtained with respect to the five measurement axes. The measurement axis of y corresponding to the minimum error is referred to as axis 0 in the next iteration, and the other new measurement axes around axis 0 are also generated using the above-mentioned method. The shrink angle α on the next iteration should be less than that on the previous iteration. . α will decrease gradually in subsequent iterations according to the shrink ratio , , so the measurement axes will close up to each other. When the shrink ratio is unchanged, the initial shrink angle determines the search scope and search speed. The initial shrink angle and the shrink ratio are specified based on the actual situation in order have a balance between the search scope and the search speed.
The method described above is based on the assumption that the principal axes of the two point clouds are fitted properly. If the principal axes are not fitted properly, the initial shrink angle must be enlarged enough to expand the search scope and the shrink ratio which is needed to be smaller so as to increase the number of iterations. Even though the angle between the two principal axes is 180° apart after the axes fitment, a good matching still will be obtained ultimately.
There must be only one center axis, even though the number of side axes is unfixed, but there must be at least three side axes.
2.3. Low-Precision Search
The measurement axes tend to limit the possible position of the x point cloud in each iteration (Figure 5). Low-precision search is used to seek for the best orientation of the x point cloud with respect to the y point cloud along every measurement axis. The low-precision search process is shown in Figure 6. In the initial iteration, the principal axis of the x point cloud coincides with axis 0 of the y point cloud. A 360° dial is sketched out on the normal plane of axis 0. When the x point cloud rotates through a step angle ( is the step angle of the low-precision search; here, = 1°) around the measurement axis, a distance error between the x and y point clouds is computed. When the x point cloud has rotated by 360° around the measurement axis, there will be 360 distance errors, and the minimum error among them will correspond to the optimal matching position of the x and y point clouds along the measurement axis (like Figure 7).
2.4. High-Precision Search
In order to obtain a better search result, a high-precision search is carried out based on the low-precision search. The low-precision search result provides an angular position with respect to the minimum error on some measurement axis. Then, in the ±n° range of this angle position (in this paper, = 1°), the x point cloud is rotated around the current measurement axis with a smaller step, Stepac (in this paper, Stepac = 0.04°), and the distance error between the x and y point clouds is calculated. Two errors per Stepac will be generated, between which the minimum is the result of the high-precision search (Figure 8).
2.5. Obtaining the Optimal Result and Subsequent Iteration Steps
The goal of the low- and high-precision search is to sort out the best positions of the x point cloud corresponding to the y point cloud along all measurement axes, as shown in Figure 9. Firstly, the principal axis of the x point cloud coincides with axis 0 of the y point cloud in accordance with the dial sketched out on the normal plane of axis 0. Then the x point cloud is rotated around axis 0, and the low- and high-precision searches are executed in sequence. Secondly, the principal axis of the x point cloud is made to coincide with axis 1 of the y point cloud, with a dial being drawn on the normal plane of axis 1, and the low- and high-precision searches are repeated until the five axes have been processed. In this approach, five measurement axes are created during every iteration, so that after the x point cloud has finished rotating around all the measurement axes, there will be five better positions of the x point cloud corresponding to the five measurement axes. The optimal solution of the current iteration corresponds to the minimum error among the five better positions (Figure 10).
The search scope will be reduced by closing up of the measurement axes in the next iteration. Here, the measurement axes required in the second iteration are obtained with shrink ratio . In the second iteration, axis 0 coincides with the measurement axis corresponding to the best position of the x point cloud from the first iteration. The optimal positions of the x point cloud corresponding to axis 1 in the first iteration and to the measurement axes in the second iteration are shown in Figure 11. Axis 0 in the second iteration coincides with axis 1 in the first iteration, and the shrink angle in the second iteration is smaller than the shrink angle in the first iteration. The search accuracy is much improved in a smaller search scope. The measurement axes in the third iteration can be determined in a similar fashion as shown in Figure 12. Assuming that the optimal solutions from the first, second, and third iterations correspond to axis 1, axis 2, and axis 1, respectively, the farthest point of the x point cloud will move on a spherical surface during the iterations. If the spherical surface was developed into a 3D plane, the motion track of the farthest point of the x point cloud would be like that shown in Figure 13.
To add further clarification to Figure 13, the shrink angle is defined aswhere denotes the number of iterations; is the initial shrink angle; is the shrink ratio of the shrink angle, .
Equation (1) makes perfectly clear that the shrink angle between the side axis and the center axis becomes smaller with the increasing number of iterations. In Figure 13, the dots with similar patterns are gradually gathering together, and because of this phenomenon, the search accuracy is improving eventually. When , the shrink angle decreasing rate becomes smaller, while the search scope becomes larger, and more areas are searched repeatedly; however, when , the search scope becomes smaller, and convergence becomes faster, but there may be blind (unsearched) areas.
As the number of iterations becomes infinitely large, the search scope is then determined by the initial shrink angle and the shrink ratio. Using the method that represents the central angle by the planar distance (Figure 13), the largest search scope is given by
When , then (2) changes to
Equation (3) can be used to determine the relationship between the initial shrink angle and the shrink ratio.
For example, assume that , °; then = 60°. This means that under these conditions, the proposed method can deal with a scenario in which the angle error between the principal axes of the x and y point clouds is less than 60°. Assuming that = 90°, that is, the initial four side axes are perpendicular to the initial center axis, this is a worst case scenario, where the principal axes of the two point clouds are directly opposite to each other, that is, = 180°. The proposed method can handle such worst case scenarios satisfactorily without searching for dead angles.
3. Experimentation and Results
Experiments have been performed with both simulated and measured point clouds to test and verify the method proposed in this paper. The simulated point clouds were used to describe the search process as well as to see the effect of initial shrink angle on the final search results. Measured point clouds from different LiDARs were used to test and verify the robustness of the method.
3.1. Simulated Point Cloud Test
The x and y point clouds including four points with their initial states are shown in Figure 14. According to the method described above, the two point clouds are aligned together when 10 iterations of geometric fitting, principal axis fitting, low-precision search, high-precision search, and optimal result computation have been carried out. Every iteration yields five errors between the x and y point clouds (the x point cloud relates the five measurement axes of the y point cloud). After completion of the 10th iteration, the 50 error results are obtained which can be used to describe the convergence of the measurement axes and the search process. The 50 errors with ° are listed in Table 1.
As shown in Table 1, five errors values, which are described as the distances between the x and y point clouds after rotating about the current measurement axis, are produced in every iteration. After the first iteration, the absolute values of errors for axes 0–4 are 1.254601, 4.110507, 5.029510, 4.575655, and 4.08057, among which the error value 1.254601 for axis 0 is at the minimum level, and therefore axis 0 corresponds to the optimum solution. It means that the center axis 0 of the five measurement axes in the second iteration must coincide with the optimal solution, that is, axis 0 taken from the first iteration. Apparently, the optimal values of the initial six iterations all correspond to center axis 0 and the same pattern is unchanged until the 7th iteration, where a smaller error is generated on axis 4. After setting the initial shrink angle to 30° (large enough) and the shrink ratio to 1/1.7, it is observed that the principal axes of the point clouds are found to be chosen quite accurately. However, this initial shrink angle and shrink ratio are relatively large leading towards a larger search scope and lower convergence speed. If the shrink angle between initial axes is reduced to 15° while keeping the other parameters unchanged, the calculated results are shown in Table 2.
Table 2 shows that the most accurate result occurred in the 5th iteration. The final result after the 10th iteration was found out to be 1.234192, which was more accurate than the previous error result of 1.234646 as shown in Figure 15. A comparison between the two results incorporating different shrink angles is drawn in Figure 16. When the initial shrink angle was 15°, the error converged more quickly.
The two tests are done using the simulated point cloud shown in the Figure 14, one is using the principal axis fitting, the other is not using the principal axis fitting. For the test’s demands, the initial shrink angle between measurement axes was set a little bigger value in order to extend the search range. According to the above discussions, while the shrink rate , the initial shrink angle °; there is no search blind angle, so the parameters are used in this experiment. After ten iterations, the errors of the algorithm with the principal axis fitting are shown in Table 3, and then the errors of the algorithm without the principal axis fitting are shown in Table 4 and the search matching accuracy comparison and analysis chart with the principal axis fitting and without the principal axis fitting is shown in the Figure 17.
From Figure 17 seen, while including the principal axis fitting, the little errors are corrected in the search matching process on the measurement axes. Without the principal axis fitting, the large errors are corrected by the measurement axes matching. Therefore, the principal axis fitting is to reduce computation cost and improve the search matching accuracy.
3.2. Measured Point Cloud Registration Test
In order to verify the proposed algorithm, the measured point clouds from MStar8000 [18] (maximum range 80 m) and Riegl LMS-Z420i (maximum range 1000 m) were used in the testing procedure.
Here, we manually extract the key points which cover the same region, and the key points are some corresponding points of the two point clouds. Structural feature points (i.e., corner points and edge points) are relatively stable feature points which should be selected as the key points. These key points are used to establish a MASM and align the two clouds with an overlap.
3.2.1. Registration Test for Two Point Clouds with a Small Overlap
In the first step, the data of the car head captured by MStar8000 were chosen for the test (in Figure 18). Two measuring sites were established from different directions and different angles so that the two point clouds were very different involving a small overlap. Five pairs of key points were selected from the two point clouds, for which absolute correspondence was not necessary, only general correspondence was fairly enough. Here, only five pairs of key points attend to computing. After 10 iterations, convergence was achieved quickly with a satisfactory registration result, as shown in Figure 19. Although there was a small overlapped region in the two point clouds, the registration result is quite satisfactory, without the “sticking-up tail” away from the common region.
There are some yellow points in Figure 20; they show the vertices of the measurement axes. Five measurement axes are created in each iteration; these measurement axes have gradually drawn close during iterations.
3.2.2. Registration Experiment for Two Point Clouds with Different Resolutions
RieglZ420 was used to capture LiDAR data for an outdoor scene. Because of the long distances between adjacent sites, the two point clouds have different resolutions in the overlap (Figure 21(a)). An ICP algorithm was given a try to align the two point clouds, but a good match was not achieved due to the gap as shown in Figure 21(b). This phenomenon occurred because some false corresponding points led the algorithm to a local minimum and ICP need a good initial coarse alignment of two point clouds. Ten pairs of key points were selected from the two point clouds to establish a MASM using the proposed method, five measurement axes are created in each iteration. MASM can avoid local minimum and align the two point clouds well without any gaps after 15 iterations (Figure 21(c)). Figure 21(d) shows the registration result of four point clouds using the proposed method. The goal of MASM is to search a better registration position between the two point clouds; the results of MASM are not affected by different resolution of the point clouds.
(a)
(b)
(c)
(d)
4. Conclusions
The proposed MASM has used sparse clouds of 3D key points instead of the original raw scans to register them. During the registration iterations, the position of the point cloud is continuously adjusted by principal axis fitting and measurement axis searching, and the point cloud can be quickly confined to the best registration position. Principal axis fitting improves computational performance, while multimeasurement axis search avoids wrong registration. Experimental results on multiple test data sets have shown global geometrical consistency of the registered point clouds.
If a more accurate registration is demanded, once MASM has been solved, Lidar scans can be fine registered with the traditional methods like ICP.
Future research is expected to continue in upgrading the manual extraction of the key points to the automated one. Further study should put focus on the effects of various kinds of the key points on the registration results to select a proper type.
Competing Interests
The authors declare no conflict of interests.
Authors’ Contributions
Shaoxing Hu and Aiwu Zhang conceived and designed the research; Sicheng Liu performed the research; Aiwu Zhang wrote the paper.
Acknowledgments
This work is supported by program for National Natural Science Foundation of China (NSFC 41071255). This work is funded by the Importation and Development of High-Caliber Talents Project of Beijing Municipal Institutions (CIT&TCD20150323) and Specialized Research Fund for the Doctoral Program of Higher Education (SRFDP20131108110005).