TRAJECTORY-BASED REGISTRATION OF 3 D LIDAR POINT CLOUDS ACQUIRED WITH A MOBILE MAPPING SYSTEM

Thanks to a hybrid georeferencing unit coupling GNSS and IMU sensors, mobile mapping systems (MMS) with lidar sensors provide accurate 3D point clouds of the acquired areas, mainly urban cities. When dealing with several acquisitions of the same area with the same device, differences in the range of several tens of centimeters can be observed. Such degradation of the georeferencing accuracies are due to two main reasons: inertial drift and losses of GNSS signals in urban corridors. The purpose of this paper is therefore to correct these differences with an accurate ICP-based registration algorithm, and then to correct the MMS trajectory using these retrieved local transformation parameters.The trajectory loop information plays a key role for that purpose. We propose a four-step method starting from a 3D point cloud with overlapping parts, and the trajectory of the mobile mapping system. First, a polygonal approximation of the trajectory is computed in order to first divide the whole registration problem in local sub-issues. Secondly, we aim to find all the potential overlapping acquired areas between these segments using simple bounding box intersections. Thirdly, for each pair of overlapping areas, an efficient variant of the ICP algorithm is proposed to (1) prune cases where segments do not share point clouds of the same areas and (2) retrieve the transformation parameters, for real overlapping cases. Finally, all these transformations are linked together, and fed into a global distance compensation problem, allowing to adjust the MMS trajectories for several passages. As a conclusion, this method is successfully applied to data acquired over Paris (France) with the Stereopolis mobile mapping system.


INTRODUCTION
Cartography has been evolving for several decades following informatics tools developments.Nowadays, 3D models have become common place in diverse application areas such as navigation, tourism, heritage, infrastructure monitoring or urban planing.Many Mobile Mapping Systems (MMS) were developed to provide accurate 3D point clouds and high-resolution images mainly from urban cities.Those systems allow the creation and the update of very fine 3D models by coupling lidar system and HD digital cameras.Moreover, aerial data are simultaneously used in order to refine 3D models reconstruction (Frueh and Zakhor, 2004) and provide more realistic models.Mobile Mapping Systems use hybrid georeferencing unit coupling GNSS, odometer and IMU (Inertial Measurement Unit) sensors in order to determine a precise position of the lidar sensors and cameras (Grinstead et al., 2006).Due to urban corridors, GNSS signals losses are quite frequent, degrading the quality of the compensated trajectory.So that, with several acquisitions of the same area, differences in the range of several tens centimeters can be observed, which does not allow accurate applications as metrology or change detection.An analysis of such errors is proposed in (Ellum and El-Sheimy, 2002).Nevertheless, this problem is not as complex as Multi-View Stereo (Frahm et al., 2010), where most of the time no georeferencing data are available.Since, overlapping areas may exist within datasets, (Clemente et al., 2007) introduce loop information in order to improve the MMS trajectory estimation.Such a redundancy allows to compute accurate registration between two point clouds of the same overlapping area.Moreover, the ICP algorithm provides a robust and accurate registration of two coarse-aligned point clouds (Chen and Medioni, 1992).Many variants of this algorithm have been propose (Rusinkiewicz and Levoy, 2001).This registration step can be for instance improved using geometrical features proposed in (Demantké et al., 2011).
Therefore, the aim of this paper is to propose a correction strategy of georeferencing inaccuracies that may appear after the acquisition process with Mobile Mapping System.The method is purely based on the analysis of 3D point cloud acquired with a single lidar device.Trajectory loop information and overlapping areas are used to compute local transformations, which are fed into a global distance compensation problem, thus allowing to adjust the MMS trajectory.After introducing the MMS trajectory estimation, the four-step method is developed in Section 3 as follows : polygonal approximation of the trajectory, research of overlapping areas, registration of each of these areas and global minimization of the trajectory.The dataset is presented in Section 4. Finally, results of the method are presented in Section 5, and conclusions are drawn in Section 6.

Estimation
Various on-the-fly systems, like GPSVan from Center of Mapping at the Ohio State University (Novak, 1991), GEOMOBIL from ICC (Alamús et al., 2005), Cyclomedia (Swart et al., 2011), Stereopolis from IGN (Paparoditis, 2011) have been developed mainly for street mapping applications.In these systems, the trajectory is estimated thanks to the direct georeferencing technology.Satellite positioning and inertial navigation information are provided by GNSS and IMU, respectively.Both technologies are integrated so that position is given by the GNSS receiver, and orientation by the IMU device.Such an integration allows to overcome their individual limitations.Since IMU provide attitude data with a higher rate than GNSS positions, IMU is used to detect and correct GNSS cycle slips and signal gaps.Conversely, IMU accelerometer calibration can be performed thanks to the GNSS.
An odometer is traditionally added to provide vehicle body velocities and minimize the errors during GNSS signal losses.Finally, the trajectory is optimally estimated with a Kalman Filter, that compensates the GNSS/IMU system errors (Abuhadrous et al., 2003).Several solutions, based on more or less complicated mathematical models exist, depending on the quality of the IMU.For more details, please see (El Sheimy, 2000).However, the IMU drift and the GNSS signal losses are such that differences are still visible after the compensation step (Figure 1).

Improvement
Many solutions have been proposed to improve the quality of the estimated MMS trajectory.As mentioned above, the first solution consists in improving the mathematical model introduced in the Kalman Filter (KF) solution.Such strategy assumes that both typical and large errors are well known, which is, in practice, not straightforward to have.The KF solution can be improved with the introduction of additional sensors or constraints.For example, new IMUs or GNSS antennas can be added, and coupled in a cascaded strategy (Cazzaniga et al., 2007).A coarse 3D model is introduced in (Lothe et al., 2009) for real-time applications.Finally, other methods directly rely on the acquired data.Consequently, both photogrammetric and lidar-based solutions are possible.Firstly, data resulting from Structure-from-Motion and dense matching are often used.Tardif et al. (2008) have noted that large field-of view cameras should be preferred to reduce the trajectory drift.Secondly, unlike photogrammetric techniques, lidar solutions are often not based on redundant measurements.Calibration is for example performed using planar patches detected of the 3D point clouds.Howard et al. (2004) propose a threestep approach in order to compute trajectory.Firstly, a "Fine-Scale" localization step use IMU and a "scan matching" point cloud registration method in order to compute accurate segment of trajectory.Then, they perform a "Coarse-Scale" localization step of each previously computed segment using GNSS.Finally, "Coarse-to-Fine" localization aims to compensate the whole trajectory.Feature matching are introduced in this last step.

Proposed strategy
Overlapping point clouds can indeed to locally registered to reduce the existing shifts.The retrieved transformations allow to introduce new constraints in the estimation procedure, and therefore may benefit to areas that have been only acquired once.The four steps are as follows: • Polygonal approximation: The trajectory is fragmented in segments in order to subsequently perform piecewise registration; • Overlapping area detection: Segments corresponding to the same geographical areas are matched; • Registration: Pairs of segments (i.e., point clouds) are registered.An improved ICP method is proposed for that purpose.
• Global trajectory estimation: Local transformations are fed into a global minimization procedure.

Polygonal approximation
In order to fragment the initial trajectory in a reduce number of segments, the method described in (Lowe, 1987) is used.A set of points is recursively cut into segments with respect to both their length and deviation.Here, the points correspond to the reconstructed trajectory i.e., the estimated MMS positions.Since the KF procedure is mastered by the georeferencing device with the highest acquisition rate, one position is available for each IMU measurement (100 Hz).This iterative algorithm consists of the following steps, illustrated in Figure 2.For a set of points {P1, ..., PN }, the segment l = P1PN is firstly considered.The farthest trajectory point from this segment is considered (Figure 2a).(2) Furthermore, the same method can be computed on the subset {P1, ..., Pi m } and {Pi m , ..., PN } (Figure 2c).This allows to define (Ψ0, Ψ1, Ψ2) as follows: In order to determine if the sub-segmentation is better than the initial one, one can check whether: where α is an adjustment coefficient that defines the level of generalization.If this inequality is true, point Pi m is conserved, and the process is recursively applied on the two sub-segments.Otherwise, segment l is conserved.Finally, the polygonal approximation can be defined by a subset of points of the trajectory {Pa 1 , ..., Pa M } where a1 = 1 and aM = N .
Similarly and subsequently to the trajectory, the lidar point cloud can be segmented.Since each 3D point is associated to a MMS position during the KF procedure, the partition is straightforward.

Overlapping area detection
In order to register areas acquired several times, overlapping parts have to be found within the dataset.A first coarse approximation can be provided by the use of axis-aligned minimum bounding box.The result is then improved using point-to-point distance.
Approximation Firstly, the bounding box Bn of each lidar subset Sn is computed.Then, for a given subset Si, non-empty intersections with all previous bounding boxes of the trajectory (Sj, i < j − ) are computed, resulting in the set Ω: The parameter allows to avoid pairs of consecutive segments.It is simply set to 2 in our method.
The following step then consists in merging segment pairs as (Si 1 , Sj), (Si 2 , Sj), when Si 1 and Si 2 are consecutive segments of the trajectory: A new set of pairs Ω f is therefore obtained.Such procedure allows to avoid keeping several consecutive small segments in order to perform robust registration.
Refinement For a given pair of lidar point sets (A, B) ∈ Ω f , matching points between A and B are retrieved.Such matching points are retrieved using simple nearest neighbor query.The number of corresponding matching points with a relatively small distance (around the point cloud resolution) is noted N (A, B).Finally, only pairs (A, B) where N (A, B) > T are conserved, resulting in the set Ωr. T is set to 100 in order to achieve robust results in the following registration step.

Registration
This step consists in registering each pair of lidar point sets.In the literature, the most commonly adopted solution is the well-know ICP (Iterative Closest Point) algorithm (Chen and Medioni, 1992;Besl and McKay, 1992).This is a robust method that has been incremented with many variants (Rusinkiewicz and Levoy, 2001).This standard approach performs a fine registration of two overlapping point clouds by iteratively estimating the transformation parameters, assuming that a good a priori alignment is provided.In case of datasets acquired with a MMS equipped with GNSS and IMU, such condition is fulfilled.Furthermore, this algorithm can be easily improved and speed up using features of interest (Sharp et al., 2002) and (Bae and Lichti, 2008).We follow this strategy, and introduce in the ICP procedure the geometrical features proposed in (Gressin et al., 2012).
Features of interest Many neighborhood descriptors exist in the literature (Sharp et al., 2002) and (Rabbani et al., 2007).Demantké et al. (2011) have proposed a multi-scale method that use geometric features computed from the local covariance matrix in order to retrieve the optimal radius neighborhood of each point.For a given radius r and his associate neighborhood V r , a principal component analysis give three eigenvalues (λ1, λ2, λ3).The standard deviation along an eigenvector i is denoted by σi = λj.Then, three geometrical features are introduced in order to describe the linear (a1D), planar (a2D) or scattered (a3D) behaviors within V r : The dimensionality labeling (1D, 2D, or 3D) of V r can be defined by: If σ1 σ2, σ3 , then a1D is greater than the two others so that the dimensionality labeling d * (V r ) results to 1. Conversely, if σ1, σ2 σ3, a2D prevails.Finally, σ1 σ2 σ3 implies d * (V r ) = 3.One example of dimensionality labeling is presented in Figure 4. Dimensionality features are computed for increasing radius values between a lower bound rmin to an upper bound rmax, using a square factor.For each radius r and for each point P , a measure of unpredictability is given by the Shannon entropy of the discrete probability distribution (a1D, a2D, a3D): Then, the optimal neighborhood radius r * is selected so as to minimize the entropy feature (cf. Figure 5).
The optimal neighborhood V * , associated to r * is then used to compute the dimensionality labeling d * (V r * P ).The more E f , the more prominent the retrieved local geometrical behaviour.

Local registration
The four steps of the ICP algorithm are: (1) selecting a reduced number of points, (2) finding matching points, (3) weighting each corresponding pair, and, finally, (4) designing and minimizing a given error metric.Using geometric features can improve several steps of this method.According to (Gressin et al., 2012), the selection and the rejection steps can be enhanced as follows.Firstly, a selection of points with a high entropy value E f (more than 0.7) can better the convergence speed of a factor of two while maintaining an accuracy of one tenth of the point cloud resolution.This corresponds to points with a clear local prominent behaviour, i.e. in urban areas, building facades (cf. Figure 5).Then, rejecting 70% of the farthest matched points using a specific distance improves the registration accuracy to a few millimeter (depending on the point cloud resolution).The following distance dV (P1, P2) = VP 1 − VP 2 has been introduced, where V = i σi is proportional to the ellipsoid volume, and allows to characterize the shape of the 3D point neighborhood.Finally, this algorithm is applied to each pair of subset point clouds (Si, Sj) ∈ Ωr.This step provides two information for each couple: 1.A transformation (T ICP S i ,S j ), which allows to register the two overlapping point clouds; 2. An indicator of the quality of the procedure (σ ICP S i ,S j ), providing by the average distance of matched points after registration.
The first one will be used for the next stage of global minimization.

Global compensation
The two previous steps linked overlapping areas.The aim is now to compute the trajectory with more constrains in order to improve its accuracy as well as the 3D point cloud accuracy.Because the positional system data (IMU, GNSS, odometer) are some times not available, our method only necessitates the 3D georeferencing position of the lidar scanner in order to obtain two different types of constraints (absolute and relative positioning constraints).Since the rotation resulting of the registration is very low, only translation-based minimization has been formulated and applied.Furthermore, introducing rotation in the minimization step results in loosing the linearity of the problem.
Let {Pa 1 , ..., Pa M } be an approximation of the trajectory {P1, ..., PM } of the mobile mapping system.Let {T ICP a i ,a j } be a set of transformation given by the previous steps.The equations that the solution {P * a i } must satisfy are : • Absolute positioning constraint : • Relative positioning constraint : • Registration constraint : This over constrained and linear problem can be solved using a standard least squares method.Furthermore, different weighting can be applied for the different equations.Indeed, the accuracy of the various positioning devices have different orders of magnitude.Whereas GNSS have typically an absolute accuracy of a few meter, the IMU and the odometer can give a relative accuracy of a few centimeter, and the ICP gives a relative accuracy of a few millimeter.So weights ranging in the same proportions are introduced.
Finally, the transformation retrieves for each point of the approximate trajectory is linearly distributed on all the points of the trajectory, and transferred to the corresponding lidar points.

DATASET
Our method have been tested on one urban dataset, over Paris (France), acquired with the Stereopolis (IGN) mobile mapping system.This dataset covers a surface of 0.5 km 2 , with a trajectory of 1.5 km along several streets.The whole point cloud includes 5 millions of points.The acquisition procedure of this dataset is based on a RIEGL LMSQ-120 lidar sensor which provides vertical profiles perpendicular to the platform trajectory.The distance range of the system is from 2 to 150 meters with an accuracy of 20 mm.The lidar acquisition frequency is up to 30 kHz.The RIEGL sensor allows a vertical coverage of 40 deg with an angular resolution of 0.04 deg, using a rotating polygon mirror.Positioned vertically, the lidar sensor allows to acquire both the roads, the sidewalks and medium-sized buildings.However, the tops of highest building may not be visible, depending on the distance from the sensor to the building.This configuration allows to cover areas with several interesting objects (street furniture), but also included many moving objects as pedestrian or cars.Using geometrical features permits to discard such items from the ICP procedure.The positional system coupling two GPS and one Applanix IMU POSVL200.The GPS acquire position at 1Hz, when the IMU have a frequency rate of 100Hz.This IMU has a drift of 3 deg/hr, resulting an accuracy of about 20 cm for a GPS signals looses of 1 minute (2 m for 3 minutes).

RESULTS
An upside view of this area can be seen in Figure 6 (left).Figure 6 (right) shows the polygonal approximation of the trajectory.This polyline is composed of about 550 segments.This number of segment, directly depend on the choice of the α parameter (Section 3.2).In this case, the value α = 0.5 is chosen in order to have enough segments, allowing a good approximation of the trajectory.This value also permits to not have too small segment that not allow to perform accurate registration.After the refinement procedure, the second step provides 20 potential overlapping areas.All these selected areas have then successfully registrated using our modified ICP algorithm.The global compensation method delivers an adjusted trajectory.This correction ranges from -30 cm to 50 cm for the vertical component (Figure 7).The planimetric correction ranges from 0 to 45 cm.Peaks visible on the graph (1, 2, 4, 6) correspond to the overlapping areas, where ICP constrains are applied.Contrariwise, not compensated areas (3, 5) correspond to the area without redundancy.Finally, this correction is propagated to the 3D lidar point cloud (Figure 8). Figure 9 illustrates the results of the procedure for three areas of interests.In these areas, significant elevation differences are visible on the input data.These degradations are correctly rectified by our method.One can see that the trajectorybased correction of the georeferencing process is satisfactory and is not limited to areas with multiple surveys.

CONCLUSION
The objective of this paper was to propose a method for correcting georeferencing inaccuracies that may appear after the lidar acquisition process with a Mobile Mapping System in a urban environment.The strategy was based on the trajectory loop information and used areas acquired several times to find the local transformations, afterwards applied to the whole area.The only mandatory input is the trajectory estimated with a standard Kalman filter.The strategy did not require as input any measurement provided by the GNSS, the IMU or the odometer.The MMS trajectory is first fragmented in order to temporally split the mapped environment into sub-areas.Each segment is examined so as to find whether overlapping areas exist.Then, the ICP algorithm is applied on point clouds belonging to pairs of segments.In order to achieve very accurate and fast registrations, the ICP is improved by the introduction of features, focusing the process on the most reliable in terms of local geometry.The obtained transformations are finally fed into a global compensation procedure that allows to diffuse the retrieved information to the other parts of the trajectory.The proposed method is successful despite the single use of translation in the compensation step.This is due to the fact that georeferencing errors are not very important in rotation.More accurate trajectories and point clouds are now available, offering the opportunity of processing point clouds with higher densities and quality or opening the field of change and mobile object detection.As the ICP algorithm can also provide a rotation term, improvement may focus on introducing an orientation term in the global compensation step.Moreover, since most of Mobile Mapping Systems also include digital cameras with HD (panoramic) images, an approach integrating both photogrammetric and lidarbased constraints is conceivable.For instance, Cannelle et al. (2009) proposed an image-based method that allows to adjust the MMS trajectory.These position data could be merged with those from our method in order to improve the reliability of the final compensation.Eventually, another possibility is to directly use, in the compensation step, raw data from the georeferencing unit that contains, as for the registration step, uncertainty information.

Figure 1 :
Figure 1: Georeferencing errors between two point clouds (red and blue color).The two trajectories are also visible.Shift can be noticed on cars and buildings facades.

Figure 3 :
Figure 3: Result of polygonal segmentation and the associated partition of the lidar point cloud over an area of interest (one color per segment, top view).One can note that all the segments do not have the same length.

Figure 6 :
Figure 6: Trajectory fragmentation.Left: raw data, colored with respect to the distance from the MMS (low (black)→high (white)).Right: extracted polyline, colored with respect to the acquisition time (blue→red).

Figure 7 :Figure 8 :
Figure 7: Results of the global compensation procedure.Areas of interest are indicated in Figure 6

Figure 9 :
Figure 9: Results of our approach for three areas of interest.Left: before registration.Right: after registration.