ACCURATE REGISTRATON OF MMS POINT CLOUDS OF URBAN AREAS USING TRAJECTORY

Recently, by Mobile Mapping System (MMS) with laser scanners, a GPS and IMU (Inertial Measurement Unit), 3D point clouds of urban areas (MMS point clouds) are easily acquired. When the same areas are scanned several times by the MMS, the point clouds often have differences in the range of several hundreds of millimetres. Such differences are caused by inertial drifts of IMU and losses of GPS signals in urban areas. In this paper, we propose an automatic accurate registration method of MMS point clouds using a new variant of ICP (Iterative Closest Point) algorithm for MMS point clouds and trajectory modification. Our method consists of four steps. Firstly, some trajectory points are automatically extracted by analyzing the trajectory. Secondly, the differences of point clouds are derived at the extracted trajectory points in the overlapping scan region by our new ICP algorithm which minimizes pointto-plane and point-to-point distances simultaneously and filters incorrect correspondences based on a point classification by PCA (Principle Component Analysis). Thirdly, the modified positions and rotation parameters at all extracted trajectory points are derived by a least squares method for positioning and registration constraints. Finally, each point in the point clouds is modified by coordinate transformations which are derived from linear interpolation of the modified positions and rotation parameters of the extracted trajectory points. Our method was applied to MMS point clouds and trajectory and the performances were evaluated.


INTRODUCTION
Mobile Mapping System (MMS) with laser scanners, GPS and IMU (Inertial Measurement Unit) contributes easy acquisitions of 3D point clouds of urban areas (MMS point clouds).When the same areas are scanned by MMS several times, MMS point clouds often have differences in the range of several hundreds of millimetres.Such differences are caused by inertial drifts in IMU and losses of GPS signals in urban areas.Therefore, an automatic registration method which accurately and effectively modifies these differences is strongly required.
Many point cloud registration algorithms have been proposed, for example 4PCS (Aiger et al., 2008) and ICP (Iterative Closest Point) algorithm (Besl andMcKay, 1992, Chen andMedioni., 1992).The ICP is especially effective for correcting the differences in MMS point clouds, because magnitude of the differences is relatively small compared with the size of point clouds and the point clouds have similar poses.ICP algorithm provides an accurate registration between two point clouds by iteratively minimizing registration errors which are squared distances between corresponding points in each point clouds.Many variants of ICP about the selection of corresponding points and error metrics have been proposed (Al-Durgham et al., 2011).Rusinkiewicz et al. (Rusinkiewicz and Levoy, 2001) evaluated some variants of ICP algorithms, and concluded that ICP using point-to-plane distances for error metric (point-toplane ICP) is more accurate than the standard ICP which minimizes point-to-point distances.Ridene et al. (Ridene and Goulette, 2009) propose a registration method for MMS point clouds.In their method, MMS point clouds are divided into blocks and point-to-plane ICP is applied to each block pair.However in this method, gaps between neighboring blocks may occur.In addition, when their method is applied to multipath point clouds, accumulated errors sometimes occur.Gressin et al. (Gressin et al., 2012) propose a MMS point cloud registration method based on trajectory.In their method, firstly, the differences are derived by an ICP algorithm at several trajectory points.Next, MMS trajectory is modified by the least squares method based on position and registration constraints.Finally, point clouds are registered according to the modified trajectory.Therefore, this registration method provides a better result because it performs continuous registration of point clouds along the trajectory.However in this method, point-to-point ICP algorithm is used, and the rotation is not considered in the registration process.Moreover, the accuracy of their method has not been evaluated.
In this paper, we propose an automatic accurate registration method for MMS point clouds using a new variant of ICP algorithm for MMS point clouds and trajectory modification and we evaluate the registration accuracy of our method.Our method follows three extensions from the method proposed by Gressin et al. (Gressin et al., 2012)

Overview of the proposed registration method
Figure 1 shows the overview of our method.Our method consists of the following four steps.
Step 1 : Difference Modification Points extraction Difference Modification Points (DMPs) for trajectory modification are automatically extracted from all the given trajectory points  by analyzing accelerations and angular velocity of MMS and trajectory intersections.This step reduces the number of redundant difference calculations in the step2.
Step 2 : Difference calculation by CCICP algorithm The differences at each DMP in the overlapping regions are derived by the CCICP (Classification and Combined ICP) algorithm.The CCICP algorithm is a new variant of ICP algorithm for MMS point clouds which minimizes point-toplane and point-to-point distances, simultaneously, and also filters incorrect correspondences based on point classification by PCA (Principle Component Analysis).
Step 3 : Derivation of transformation parameters Translation and rotation parameters for all DMPs are derived as least squares solutions for absolute positioning constraint, relative positioning constraint, and registration constraint.
Step 4 : Registration of point clouds Each point position   in the point cloud is modified by coordinate transformations which are derived from linear interpolation of the modified positions of DMPs and the derived rotation parameters.

Difference Modification Points extraction
The differences often occur at areas where MMS velocity and orientation change drastically.Therefore, in these areas, the differences are non-linear and point clouds may be distorted.In addition, point clouds at intersecting trajectory regions should be registered accurately.From these points of view, in our method, DMPs are automatically extracted from all given trajectory points  by analyzing accelerations and turning speeds of MMS and trajectory intersections.DMP extraction consists of following two steps.
Step Step 2 : Difference calculation by CCICP Step 3 : Derivation of transformation parameters Step 4 : Registration of point clouds Step

Corner points
Trajectory points Candidates of DMP

CCICP algorithm
The features of the CCICP algorithm are minimizing point-toplane and point-to-point squared distances simultaneously and filtering incorrect correspondences based on point classification by PCA.These features are shown as follows.

A) Filtering incorrect correspondences
The points in the local point clouds are classified into linear points, planar points and scatter points depending on the results of the PCA.In the ICP, linear-planar and scatter-planar correspondences are rejected as incorrect correspondences for accurate registration by the CCICP algorithm.For example, corresponding points on electric cables and building facades or ones on utility poles and walls are not used in the ICP.

B) Minimizing point-to-plane and point-to-point distances
Point-to-plane and point-to-point distances are minimized simultaneously.Point-to-plane distance minimization is applied to planar-planar correspondences, because reliable normal estimations can be performed at planar points.Point-to-point distance minimization is applied to the other correspondences.Therefore points on roads, building facades and walls are registered accurately, and points on utility poles and electric cables are also registered.In this paper, point-to-plane distance minimization problem is linearized and solved using the method of Low (Low, 2004).
The CCICP algorithm follows five steps.In the following algorithm,  and T indicate a source point cloud and a target point cloud.   and    indicate a position of point  ∈  and a position of point  ∈ .
Step 1 : Selection Subset of  for distance calculation and rigid transformation derivation is selected.The selection is done by random sampling (  [%] of  are selected).The sampled point set is represented by '.
Step 2 : Point classification Each point  in ′is classified into linear points, planar points and scatter points by PCA (Demantke et al., 2011).In this method, local point distributions of each point  are evaluated by eigenvalues  1  ,  2  and  3  (  1  ≥  2  ≥  3  ) and unit eigenvectors  1  ,  2  and  3  , of the variance-covariance matrix   which is calculated from the position of point  and its neighbors (points within distance   from  ).In order to evaluate the local point distributions, the magnitude relations of the eigenvalues are investigated using the following three values: is larger than the others, the local point distribution of point  is recognized as linear.If  2  is larger than the others, it is recognized as planar.In this case, eigenvector  3  is kept as the normal vector of .In addition, if  3  is larger than the others, it recognized as scatter.
Step 3 : Matching The point ′ in T which is closest to point  ∈ ′ is extracted as a corresponding point of point  .If ′ is found, it is classified by PCA.
Step 4 : Rejecting Correspondences consisting of a linear point and a planar point or a planar point and a scatter point are rejected.A set of resulting correspondences is denoted by , and a set of planar-planar correspondences is denoted by .
Step 5 : Minimizing The average of point-to-plane squared distance  _ is derived from equation (1).
Where   ′ = an normal vector of point  ′ in a homogeneous coordinate system  = a transformation matrix in a homogeneous coordinate system The averages of point-to-point squared distances  _ are derived from equation (2).
From equations ( 4) and ( 5), the parameters of the coordinate transformation   = (,   ), which minimizes  _ and  _ simultaneously, are derived from equation ( 6).Where,  is a rotation parameter about the , ,  axis and   is translation parameter.
Where  ̂′ = [ ̂1 ⋯  ̂|∖| ] T (6) The solution   is obtained by solving  ′  = ′ using pseudo inverse matrix of  ′ .At each iteration, source point cloud  is transformed by the   .Here if the number of iterations from step3 to step5 reaches threshold   or the sum of  _ and  _ is less than threshold   or the difference between current and previous sums in the iterations is less than   , algorithm is terminated, otherwise return to step3.Therefore the amount of difference is determined by a concatenated transformation of   at each iteration.

Derivation of transformation parameters
In order to calculate the transformation parameters for all DMPs, the modified position and rotation parameters are derived by solving the overconstrained equations defined by absolute positioning constraint, relative positioning constraint, and registration constraint.The constraints are defined by equations ( 7)-( 9).
 Absolute positioning constraint:  Relative positioning constraint:  Registration constraint: Where   * = a modified position of   which is derived by a least square method   = a set of source DMPs   * = rotation angles around , ,  axis which are derived by a least square method

Registration of point clouds
Rotation and translation parameters (  ,   ) for point  in point clouds are derived by linear interpolation of modified position and rotation parameters of neighboring DMPs  − 1 and  ( −1 ≤   ≤   ).The linear interpolation is denoted by equation ( 10) and equation ( 11).
Where   =   * −   Finally, position of each   in the point clouds is modified by the coordinate transformation using   and   .
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey  5. Some parameters, for example  ℎ and  ℎ , were determined by common knowledge about running condition of MMS, on the other hand, other parameters were determined by experimentally .

Registration results and error evaluation
The results of our method are shown in Figures 6-8.On the other hand, in the right figures, the differences are reduced to about less than 100mm at roads and building facades.Figure 9 shows the differences at position C using another color range.In this figure, as a result, differences are less than 50mm at roads and building facades after registration.These results show that an accurate registration of point clouds which reduce the differences among point clouds to less than 50mm was realized by our method.Moreover, the accuracy of the CCICP algorithm was compared with a standard ICP algorithm (random sample, no rejection, point-to-point, and a quaternion based solution method).To evaluate CCICP algorithm, a point cloud that consisted of about 5.3M points and another point cloud that consisted of about 5.9M points were used.The registration parameters are the same as the ones shown in Table 5. Figure 10 shows changes of two RMS values of corresponding point distances at each iteration of the ICP.In this figure, the RMS of CCICP is 40% smaller than one of the standard ICP, and this means that our ICP realizes more accurate registration than the standard ICP.The processing times were 78.3sec for CCICP and 83.3sec for a standard ICP.It was estimated that the main difference of the processing times was caused by a difference of the solution methods for minimization problem.Coordinate transformation parameters of our ICP was derived by solving simultaneous equation using pseudo-inverse matrices.On the other hand the parameter of the Standard ICP was derived by a quaternion based method.
Finally, the effectiveness of the rotation was evaluated as shown in Figure 11 using the same point clouds and settings for accuracy evaluation.Figure 11 shows RMS values which were derived from the sum of  _ and  _ .In this figure, the RMS of CCICP with rotation is smaller than that of its CCICP without rotation.Therefore, introducing the rotation to CCICP is effective for accurate registration of MMS point clouds.

CONCLUSION
In this paper, an automatic accurate registration method of MMS point clouds using CCICP algorithm for MMS point clouds and modifying trajectory is proposed.The performance of the method is evaluated by applying it to MMS point clouds and a trajectory.Our method can reduce the differences among point clouds to less than 50mm at roads or building facades on average.Also, its computation time was about 161sec for our test data including 140M points and 2.2K trajectory points.Moreover, the accuracy of CCICP algorithm was also evaluated.
It is concluded that our ICP algorithm was more accurate than the standard ICP algorithm.Future work includes automatic parameter determination.

Figure 1 .
Figure 1.The overview of our method Figure 2. Selecting Candidates of DMP (StepA)

Figure
Figure 3. Determination of DMP and correspondences detection (StepB)

Figure
Figure 4. MMS point clouds and trajectoryTable 5. Used parameters Figure 6 shows DMPs which are extracted from trajectory points.The number of DMPs was 40, and that of DMP pairs was 19.The left figures in Figure 7 show point clouds before registration at the positions A-C which are shown in Figure 6, and right figures show point clouds after registration.In the left figures of Figure 7(a) and (c), horizontal differences are visible before registration, on the other hand, the differences are modified after registration as shown in the right figures.Similarly, in the left figure of Figure 7(b), significant elevation differences are visible.On the other hand, differences are modified as shown in the right figures.From these results, our method can perform visually accurate registration of point clouds.Our algorithm was implemented on a standard PC with Intel Core i7 3.30GHz CPU, 32GB RAM, GeForce GTX570 GPU.The computation time of registration was about 161sec.The most of computation time was the difference calculation by CCICP, and it was about 133sec.Figure 8 shows the amount of difference with color.In the figures, the amounts of difference are derived by calculations of distances between each point in the outward point cloud and its closest point in the backward point cloud.The color bar is shown in the top of Figure 8.If the closest point does not exist inside of 400 mm area from each point, points are colored in black.In left figures in Figure 8, original point clouds have about 200-400mm differences.
In this paper, MMS point clouds and MMS trajectory are denoted by and respectively.The point cloud is a set of laser scanned points  = {(  ,   )|  = (  ,   ,   ),  = 1, … , } and each point  has a position (  ,   ,   ) in a world coordinate system.Also, each point  is acquired at the GPS time   .
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey the set of DMPs  and the set of DMP correspondences , are extracted as shown in Figure 3(d).If one or more candidates or pairs of candidates exist in a block, one candidate or pair of candidates is extracted.Extracting local point clouds for ICP Our new ICP algorithm is applied to local point clouds which are extracted by DMPs and their GPS times.Firstly, two trajectory points position,    ,   are selected as the farthest back and forth points position from each DMP within   along the trajectory points.Next, the point  within the time range (   ≤   ≤    ) is extracted from the point clouds.