Algorithmic Solutions for Computing Precise Maximum Likelihood 3D Point Clouds from Mobile Laser Scanning Platforms

Abstract

Mobile laser scanning puts high requirements on the accuracy of the positioning systems and the calibration of the measurement system. We present a novel algorithmic approach for calibration with the goal of improving the measurement accuracy of mobile laser scanners. We describe a general framework for calibrating mobile sensor platforms that estimates all configuration parameters for any arrangement of positioning sensors, including odometry. In addition, we present a novel semi-rigid Simultaneous Localization and Mapping (SLAM) algorithm that corrects the vehicle position at every point in time along its trajectory, while simultaneously improving the quality and precision of the entire acquired point cloud. Using this algorithm, the temporary failure of accurate external positioning systems or the lack thereof can be compensated for. We demonstrate the capabilities of the two newly proposed algorithms on a wide variety of datasets.

1. Introduction

Laser scanning provides an efficient way to actively acquire accurate and dense 3D point clouds of object surfaces or environments. Mobile scanning is currently used for modeling in architecture and agriculture, as well as urban and regional planning. Modern systems, like the Riegl VMX-450 and the Optech Lynx Mobile Mapper, work along the same basic principle. They combine a highly accurate Global Navigation Satellite System (GNSS), a high precision inertial measurement unit (IMU) and the odometry of the vehicle to compute fully timestamped trajectories. Using a process called motion compensation, this trajectory is then used to georeference the laser range measurements that were acquired by the 2D laser scanner also mounted on the vehicle. The quality of the resulting point cloud depends on several factors:

• The calibration of the entire system, i.e., the accuracy to which the position and orientation of each individual sensor in relation to the vehicle has been determined.
• The accuracy of the external positioning sensors, i.e., the GNSS, IMU and odometry.
• The availability of GNSS, as it may suffer temporary blackouts under bridges, in tunnels and urban canyons.
• The accuracy of the laser scanner itself.

The non-availability of the GNSS in constricted spaces is an exceptionally limiting factor for mobile scanning. Without the requirement for a GNSS device, applications, such as construction or maintenance of tunnels and mines and facility management, would be open to the use of mobile scanning systems.

In this paper, we propose a novel algorithm for the calibration of a mobile scanning system that estimates the calibration parameters for all system components simultaneously without relying on additional hardware. We also deal with the problem of erroneous positioning of the vehicle in a novel fashion. We present an algorithm for computing a corrected trajectory that produces point clouds that are a more precise representation of the measured environment. To this purpose, we constructed a mobile platform with a modified approach to mobile laser scanning. Instead of a 2D laser scanner common for mobile scanning systems, we employ a 3D laser scanner as depicted in Figure 1a,c. The laser scanner obtains 3D data by rotating around its vertical axis during the scanning process (see [1]) at 6 s per rotation. This design retains the high degree of automation and speed of common mobile laser scanning systems. Furthermore, only a single laser range sensor is required to produce models with minimal data shadowing. The additional rotation with respect to the platform during the scanning process provides increased observability of positioning errors with respect to 2D mobile laser scanning systems. We evaluate the algorithms on datasets acquired by the mobile platform Irma3D and an automobile, as well as on data acquired by a commercial mobile laser scanning system (see Figure 1).

2. State of the Art

2.1. Calibration

In order for the mobile laser scanner to acquire high quality 3D models, the position and orientation of every individual sensors must be known. This paper is concerned with algorithmic calibration of these systems, i.e., algorithms to establish the parameters that best describe sensor displacements based on the sensor data itself. In this process, calibration parameters that are only roughly measured with external instruments are fine-tuned automatically.

The most basic sensor of the vehicle is its odometry, which estimates the 3D pose, i.e., the 2D position and 1D orientation, of the vehicle by counting and extrapolating the wheel rotations. Martinelli et al. present a method to calibrate the odometry readings using an augmented Kalman filter [2]. Their algorithm estimates the odometry parameters on the basis of pose information acquired by the Simultaneous Localization and Mapping (SLAM) system. As odometry is the least accurate estimator of a vehicle’s pose, mobile laser scanners are usually also equipped with more exact positioning sensors, e.g., an IMU or a GNSS device. Traditionally, these are calibrated against other positioning devices whose pose in relation to the vehicle is already known [3]. In [4], Nebot and Durrant-Whyte present a method for calibrating a low cost six degrees-of-freedom IMU on a land vehicle. Initially the vehicle is at rest to estimate the gravitational and other biases of the accelerometers and gyros of the IMU. Then, the redundancies in the measurements are exploited to estimate the calibration parameters of the IMU during a test drive.

Finally, the position and orientation of the laser measurement device also requires calibration. This is often done using a process called boresight alignment. Boresight calibration is the technique of finding the rotational parameters of the range sensor with respect to the already calibrated IMU/GNSS unit. Skaloud and Schaer describe a calibration method where an airplane equipped with a laser scanner flies several times over a residential area. Planes are extracted from the roofs in every flyover [5]. Then, the planes are matched against each other to minimize the calibration error. A similar method was developed by Rieger et al. for non-airborne kinematic laser scanning [6]. Here, the vehicle drives past the same house several times. Again, the planar surfaces of the buildings are exploited to estimate the calibration parameters of the laser scanner. The calibration parameters of the laser scanner can also be estimated when the vehicle itself is stationary. Talaya et al. present a calibration method for estimating the boresight parameters of a laser scanner by registering several scans of the environment at different positions [7]. The position and orientation of the vehicle are known at any one point, and the scans are registered using landmarks. Recently, Underwood et al. presented an approach for calibrating several range scanners to each other with no information about the pose of the vehicle [8]. The vehicle scans a previously known environment from several positions. The range data is then manually labeled, so that the ground truth for each data point is known and an error metric can be constructed. Minimizing the error yields optimal calibration parameters for the range sensors.

Our calibration system is similar to that in [8] in that multiple components are calibrated using a quality metric of the reconstructed point cloud. However, our approach also calibrates the odometry, and we require no manual selection of points or any special environment for the calibration. Instead, we employ a quality metric that is similar to the one used by [9].They calibrate a terrestrial 3D laser scanner of their own design by computing the minimum of the quality metric with respect to the internal calibration parameters. We also improve upon the metric by reducing the time complexity that is involved in its evaluation. Furthermore, we estimate the calibration parameters for multiple sensors simultaneously. This allows us to forgo a separate calibration process for every subsystem.

2.2. Simultaneous Localization and Mapping

Aside from sensor misalignment, a second source of errors are timing-related issues. On a mobile platform, all subsystems need to be synchronized to a common time frame. This can be achieved with pure hardware via triggering or with a mix of hard and software, like pulse per second (PPS) or the network time protocol [10]. Good online synchronization is not always available for all sensors. Olson developed a solution for the synchronization of clocks that can be applied after data acquisition [11].

An even more significant source of errors is the incorrect positioning of the vehicle. Solving this problem requires approaches other than classical rigid SLAM algorithms. An area that may provide a solution is the area of non-rigid registration, which is largely unexplored, except in the medical imaging community, where it is widespread, due to the need to align data with multiple modalities [12,13].

Williams et al. describe an extension of a rigid registration algorithm that includes point estimation to compensate for noisy sensor data [14]. This technically constitutes a non-rigid registration algorithm designed for low-scale high-frequency deformations. Similarly, Pitzer and Stiller provide a probabilistic mapping algorithm for 2D range scans, where also optimized point measurements are estimated [15].

Chui and Rangarajan proposed a point matching algorithm that is capable of aligning point clouds with each other [16]. These approaches are usually time expensive, due to the enlarged state space. Aiming to reduce noise, Unnikrishnan and Hebert locally fit high-order polynomials to 3D data [17]. However, large-scale deformations will not be corrected by such a local approach. Mitra et al. rely on a high number of scans that contain a slowly deforming object [18]. Instead of computing point correspondences, they estimate the deformation by the local space-time neighborhood. Consequently, higher point density is not a difficulty, but a benefit for the algorithm. On the other hand, the algorithm is incapable of dealing with large deformations and movements between individual scans.

Brown and Rusinkiewicz developed a global non-rigid registration procedure [19]. They introduced a novel variant of the Iterative Closest Point (ICP) algorithm to find very accurate correspondences. The original ICP algorithm [20] aligns laser scans acquired at different positions based on the distances between points from the two point clouds. Even though the registration requires extreme subsampling, the deformation can be successfully generalized to the entire scan. Unfortunately, this technique is not fully applicable to laser scans acquired by mobile robots [21].

Stoyanov and Lilienthal present a non-rigid optimization for a mobile laser scanning system [22]. They optimize point cloud quality by matching the beginning and the end of a single scanner rotation using ICP. The estimate of the 3D pose difference between the two points in time is then used to optimize the robot trajectory in between. A similar approach, by Bosse et al., uses a modified ICP algorithm with a custom correspondence search to optimize the pose of six discrete points in time of the trajectory of a robot during a single scanner rotation [23]. The trajectory in between is modified by distributing the errors with a cubic spline. This approach is extended to the Zebedee, a spring mounted 2D laser scanner in combination with an IMU [24]. Their algorithm sequentially registers a short window of 2D slices onto the previously estimated 3D point cloud by surface correspondences.

The approach presented in this paper optimizes the point cloud using full 6D poses and is not restricted to a single scanner rotation or time window. Instead, we improve scan quality globally in all six degrees of freedom for the entire trajectory. Finally, a rotating 3D scanner is not necessary for the algorithm to improve scan quality. This is demonstrated on state-of-the-art mobile laser scanners.

3. Calibration

Calibration is the process of estimating the parameters of a system. For a sensor system, this means that the closer these values are to the true physical quantities, the more accurate the final measurements of the entire system will be. We have designed a calibration process for mobile laser scanners that is capable of estimating all those quantities that influence the internal precision of the acquired point cloud. We assume some rough estimate is available for each parameter. The first step is to acquire a dataset with the mobile scanning system. Although it is not necessary for the geometry of the environment to be known, the system should exhibit linear motion, left and right turns and, ideally, even pitching and yawing motions. A specially designed measure of the internal precision of the point cloud is then used to find the calibration parameters via Powell’s method.

In the context of mobile laser scanners, there are several types of parameters of note. First, there is the geometrical alignment of each subsystem with respect to the vehicle. There are many frames of reference on a mobile platform. The challenge of establishing the transformation between the vehicle and the global reference system is subject to the measurements of the positioning systems and is discussed in the next chapter. For proper data acquisition, the full six degrees of freedom (DOF) pose W_s = (t_x,s, t_y,s, t_z,s, θ_x,s, θ_y,s, θ_z,s) of each sensor, s, with respect to the vehicle frame is essential. Here, t_x,s, t_y,s and t_z,s are the positional parameters specifying the translation along the x, y and z axis, whereas the θ_x,s, θ_y,s and θ_z,s parameters are Euler angles specifying the orientation around the respective coordinate axes of sensor s. Incorrect geometrical calibration leads to incorrect trajectory estimation and systematic errors in the final point cloud.

Aside from the general alignment parameters, there are internal parameters that are specific to certain sensors. These quantities directly influence the measurement accuracy of the respective sensor and could in principle be calibrated independently of the rest of the system. An example of this can be found in [25]. The laser measurements are dependent on the internal alignment of the emitter and the mirror geometry of the laser scanner, whereas the odometry is dependent on wheel circumference, w, and axis length, a, of the vehicle. In practice, for many sensors, the modification of the calibration parameters is unnecessary or infeasible, as the internal calibration is often performed by the manufacturers themselves.

Finally, systematic timing errors due to latencies can be counteracted by offset parameters o_s. We assume all sensor measurements are timestamped. Time frames are synchronized by an offset that represents the minimal inherent delay between a measurement and its reception in the system. In principle, the proposed algorithm is capable of adjusting every calibration parameter, including timing-related parameters. However, systematical synchronization errors are minor for mobile laser scanning systems and do not contribute to the quality of the final point cloud in a significant way. Therefore, we exclude timing-related parameters from the following formulation.

3.1. Algorithm Description

The principle behind the approach is to find the parameters, C, that produce the most consistent point cloud possible. The calibration parameters for all sensors, s, are concatenated to construct the calibration vector:

$$\mathbf{C}:=\hspace{0.17em}(a,\hspace{0.17em}w,\hspace{0.17em}{\mathbf{W}}_0,\hspace{0.17em}{o}_0,\hspace{0.17em}\dots ,\hspace{0.17em}{\mathbf{W}}_n,\hspace{0.17em}{o}_n)$$

(1)

In Equation (1), a and w represent odometry parameters, while W_s represent the rigid transformation of sensor s to the vehicle coordinate system. Then, the process of extracting the point cloud P = p₀, . . .,p_n consisting of n points with p_i = (x_i, y_i, z_i) from M, the entirety of measurements of every sensor, can be said to be a function f(M,C) = P. To find the optimal calibration, we must define an appropriate quality measure on P.

We employ a general quality measure that is similar to the one used by Sheehan et al. [9]. We model the points, p_i, as drawn from a probability distribution function (pdf) d(l), which represent the probability that a specific location, l, has been measured. The pdf can be approximated as:

$$d(\mathbf{l})=\frac{1}{n}\sum _j^{n}G(\mathbf{l}-{\mathbf{p}}_j,\hspace{0.17em}{\sigma }^2\mathbf{I})$$

(2)

where G(μ, σ²I) is a Gaussian with mean μ and covariance σ²I. A Gaussian distribution does not model errors, such as those that are introduced by uncertainties in the trajectory or that are dependent on the range, incidence angle or reflectance of the surface that was measured. However, it is more than sufficient to capture the notion of point cloud consistency. Calibration errors lead to surfaces appearing at multiple positions in the point cloud. The entropy of d(l) increases with these errors and decreases the more compact the point cloud is. Thus, an entropy measure on d(l) is also a quality measure for P. Sheehan et al. [9] is derived the following simplified entropy measure, which depends only on the pairwise distance of every possible pair of sample points:

$$E^{\prime }(\mathbf{P})=-\sum _i^n\sum _j^{n}G({\mathbf{p}}_i-{\mathbf{p}}_j,\hspace{0.17em}\hspace{0.17em}2{\sigma }^2\mathbf{I})$$

(3)

Sheehan et al. [9] use the Jacobian of E′ with respect to the calibration parameters of their system to apply Newton’s method for optimization. This is not possible in our case for several reasons. First, the calibration is supposed to be general, i.e., no definitive system to calibrate for is given. Second, the inclusion of parameters for the positioning systems makes the derivation of the Jacobian infeasible. This is due to the fact that in order to compute a global measurement, p_i, at time t_i the pose estimate, V̄_i, of the vehicle must be known. However, to compute W_i, all sensor measurements prior to t_i may be taken into account. Furthermore, the presence of multiple positioning sensors requires sensor fusion, thereby increasing the non-linearity and complexity of the entropy measure. In addition, we acquire a large number of sample points, usually in the order of several millions for properly calibrating an entire mobile platform. The quality measure E′(P) is infeasible for the calibration using large point clouds. One way of dealing with this problem is to reduce the number of points. We propose to uniformly subsample the entire point cloud. Furthermore, we propose to simplify the measure by using only a subsample of all pairs of points. For every point, p_i, that remains from the initial sub-sampling, we determine its closest point q_i ∈ P, such that |t_i – t_j| > δ. Here, δ is the minimal amount of time that must have elapsed for the laser scanner to have measured the same point on the surface again. Temporally close measurements are usually spatially close, as well, so they must be excluded to prevent them from dominating the quality measure. δ is easily derived from $\delta =\frac{2\pi }{v_s+v_r}$, where v_s is the rotational speed of the laser scanner and v_r, the maximal rotational speed of the vehicle. At most n, pairs of closest points are thus used in the evaluation of the error metric instead of n². We seek to find:

$$\widehat{\mathbf{C}}=\underset{\mathbf{C}}{\text{argmax}}\hspace{0.17em}E(f(M,\hspace{0.17em}\mathbf{C}))$$

(4)

where

$$E(f(M,\hspace{0.17em}\mathbf{C}))=-\sum _i^{n}G({\mathbf{p}}_i-{\mathbf{q}}_i,\hspace{0.17em}2{\sigma }^2\mathbf{I})$$

(5)

Standard minimization algorithms compute the derivative of the error function with respect to the calibration parameters. This is infeasible for the function, E(f(M, C)), as it involves complex algorithms for filtering data and fusing multiple modalities of measurements. We employ Powell’s method for optimizing E as it does not require the derivative [26]. Instead, we must merely provide an initial estimate of C, which is readily obtainable by manual estimation.

3.2. Strategies for Increasing Efficiency

To deal with the massive amount of data in a reasonable time, we employ several strategies. First, we uniformly and randomly reduce the point cloud by using only a constant number of points per volume, typically to one point per 3 cm³. An octree data structure is ideally suited for this type of subsampling. An additional advantage of the subsampling uniformity is that surfaces closer to the scanner do not unfairly influence the quality measure more than surfaces that are farther away. We also employ an octree data structure as a fast and compact representation of the entire point cloud. This allows us to compute the nearest neighbor for a point at comparable speed to state-of-the-art k-d tree search libraries [27]. The octree also compresses the point cloud, so it can be easily stored and processed.

We employ the nearest neighbor search algorithm as described in [28], which has been modified to enforce the time constraint. The modified algorithm is summarized in Algorithm 1.

Algorithm 1. FindClosest

Input: query point q, maximal allowed distance, d
1:	function FindClosest(q, d)
2:	lookup deepest node N containing bounding box of q
3:	convert q to octree coordinate i
4:	return FindClosestInNode(N, q, i, d, 0)
5:	end function
1:	function FindClosestInNode(N, q, i, d, p)	▹ p is the currently closest point
2:	compute child index c from i
3:	forj = 0 to 7 do
4:	get next closest child C = N.children[ sequence [c] [j]]
5:	ifC is outside of bounding ball then
6:	return currently closest point p
7:	else
8:	ifC is a leaf then
9:	FindClosestInLeaf(C, q, d, p)
10:	else
11:	FindClosestInNode(C, q, i, d, p)
12:	end if
13:	end if
14:	end for
15:	return currently closest point p
16:	end function
1:	function FindClosestInLeaf(C, q, d, p)
2:	T(q) is timestamp of q
3:	for all points o do
4:	T (o) is timestamp of point o
5:	if |T (q) – T (o)| > δ and |o – q| < d then
6:	p ← o
7:	d ← |o – q|
8:	end if
9:	end for
10:	return closest point p
11:	end function

Algorithm 1. FindClosest

Input: query point q, maximal allowed distance, d
1:	function FindClosest(q, d)
2:	lookup deepest node N containing bounding box of q
3:	convert q to octree coordinate i
4:	return FindClosestInNode(N, q, i, d, 0)
5:	end function
1:	function FindClosestInNode(N, q, i, d, p)	▹ p is the currently closest point
2:	compute child index c from i
3:	forj = 0 to 7 do
4:	get next closest child C = N.children[ sequence [c] [j]]
5:	ifC is outside of bounding ball then
6:	return currently closest point p
7:	else
8:	ifC is a leaf then
9:	FindClosestInLeaf(C, q, d, p)
10:	else
11:	FindClosestInNode(C, q, i, d, p)
12:	end if
13:	end if
14:	end for
15:	return currently closest point p
16:	end function
1:	function FindClosestInLeaf(C, q, d, p)
2:	T(q) is timestamp of q
3:	for all points o do
4:	T (o) is timestamp of point o
5:	if |T (q) – T (o)| > δ and |o – q| < d then
6:	p ← o
7:	d ← |o – q|
8:	end if
9:	end for
10:	return closest point p
11:	end function

Closest points are accepted only if |p_i – q_i| ≤ d_max. Although d_max is an input parameter that has to be manually set for the algorithm, it massively reduces the computation time needed to evaluate the quality measure. Furthermore, it is linked to the second input parameter, σ. Point pairs with the maximal distance, d_max, should contribute to the quality metric with a fixed ratio, r, to that of the “perfect” point pair with zero distance, i.e., we set $\sigma =\sqrt{-\frac{d_{\text{max}}}{2\hspace{0.17em}\text{ln}\hspace{0.17em}r}}$. In the experiments, the ratio has always been set to r = 0.01.

4. Semi-Rigid Optimization for SLAM

In addition to the calibration algorithm, we also developed an algorithm that improves the entire trajectory of the vehicle simultaneously. Unlike previous algorithms, e.g., by [22,23], it is not restricted to purely local improvements. We make no rigidity assumptions, except for the computation of the point correspondences. We require no explicit vehicle motion model, although such information may be incorporated at no additional cost. The algorithm requires no high-level feature computation, i.e., we require only the points themselves.

The movement of the mobile laser scanner between time t₀ and t_n creates a trajectory T = {V₀, . . ., V_n}, where V_i = (t_x,i, t_y,i, t_z,i, θ_x,i, θ_y,i, θ_z,i) is the 6 DOF pose of the vehicle at time t_i with t₀ ≤ t_i ≤ t_n. Using the trajectory of the vehicle, a 3D representation of the environment can be obtained by “unwinding” the laser measurements, M, to create the final map, P . However, sensor errors in odometry, IMU and GNSS, as well as systematic calibration errors and the accumulation of pose errors during temporary GNSS outages degrade the accuracy of the trajectory and, therefore, the point cloud quality.

The current state-of-the-art developed by [23] for improving the overall map quality of mobile mappers in the robotics community is to coarsely discretize the time. This results in a partition of the trajectory into sub-scans that are treated rigidly. Then, rigid registration algorithms, like the ICP and other solutions to the SLAM problem, are employed. Obviously, trajectory errors within a sub-scan cannot be improved in this fashion. Applying rigid pose estimation to this non-rigid problem is also problematic, since rigid transformations can only approximate the underlying ground truth. Consequently, overall map quality suffers as a result.

We employ a much finer discretization of the time, at the level of individual 2D scan slices or individual points. This results in the set of measurements M = {m₀, . . ., m_n}, where m_i = (m_x,i, m_y,i, m_z,i) is a point acquired at time t_i in the local coordinate system of V_i. In case V_i represents more than a single point, V_i is the local coordinate of the first point. All represented points are motion compensated with the best known interpolated trajectory. As modern laser scanners typically acquire 2D scan slices at a frequency of 100–200 Hz, time is discretized to 5–10 ms. Applying the pose transformation, ⊕, we derive the point p_i = V_i ⊕ m_i = R_{θx,i, θy,i, θz,i} m_i + (t_x,s, t_y,s, t_z,s)^T in the global coordinate frame and, thereby, also the map P = {p₀, . . ., p_n}. Here, R_{θx,i, θy,i, θz,i} is the rotation matrix that is defined by:

$$\left(\begin{array}{ccc}\text{cos}\hspace{0.17em}{\theta }_y\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{z,i}& -\text{cos}\hspace{0.17em}{\theta }_{y,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{z,i}& \text{sin}\hspace{0.17em}{\theta }_{y,i}\\ \text{cos}\hspace{0.17em}{\theta }_{z,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{x,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{y,i}+\text{cos}\hspace{0.17em}{\theta }_{x,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{z,i}& \text{cos}\hspace{0.17em}{\theta }_{x,i}\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{z,i}-\text{sin}\hspace{0.17em}{\theta }_{x,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{y,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{z,i}& -\text{cos}\hspace{0.17em}{\theta }_{y,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{x,i}\\ \text{sin}\hspace{0.17em}{\theta }_{x,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{z,i}-\text{cos}\hspace{0.17em}{\theta }_{x,i}\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{z,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{y,i}& \text{cos}\hspace{0.17em}{\theta }_{z,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{x,i}+\text{cos}\hspace{0.17em}{\theta }_{x,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{y,i}\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_{z,i}& \text{cos}\hspace{0.17em}{\theta }_{x,i}\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{y,i}\end{array}\right)$$

(6)

Given M and T, we find a new trajectory T′ = {V′₁, . . ., V′_n} with modified poses, so that P generated via T′ more closely resembles the real environment.

The complete semi-rigid registration algorithm proceeds as follows. Given a trajectory estimate, we compute the point cloud, P, in the global coordinate system and use nearest neighbor search to establish correspondences. Then, after computing the estimates of pose differences and their respective covariances, we optimize the trajectory, T. This process is iterated until convergence, i.e., until the change in the trajectory falls below a threshold.

To deal with the massive amount of data in reasonable time, we employ the same strategies as for the calibration. Using an octree data structure, we uniformly and randomly reduce the point cloud by using a constant number of points per volume, typically to one point per 3 cm³. Furthermore, in initial stages of the algorithm, estimates V̄_i,j are only computed for a subset of poses V₀, V_k, V_2k, . . ., with k in the order of hundreds of milliseconds. In every iteration k is decreased so that the trajectory is optimized with a finer scale.

4.1. Pose Estimation

Our algorithm incorporates pose estimates from many sources, such as odometry, IMU and GNSS. For this paper, we use the formulation of pose estimates, V̄_i,i+1, that are equivalent to pose differences:

$${\overline{\mathbf{V}}}_{i,i+1}={\overline{\mathbf{V}}}_i\ominus {\overline{\mathbf{V}}}_j$$

(7)

The operator, ⊖, is the inverse of the pose compounding operator, ⊕, such that V_j = V_i ⊕ (V_i ⊖ V_j). In addition to the default pose estimates that may also be enhanced by separating all pose sensors into their own estimates with proper covariances, we estimate differences between poses via the point cloud, P.

For each measurement, p_i, we find a closest measurement, p_j, via nearest neighbor search with |t_i – t_j| > δ, where δ is again the minimal amount of time that must have elapsed for the laser scanner to have measured the same point again. Points are stored in the global coordinate frame as defined by the estimated trajectory, T. The positional error of two poses, V_i and V_j, is then given by:

$$E_{i,j}=\sum _{k=i-N}^{i+N}{\Vert {\mathbf{V}}_i\oplus {\mathbf{m}}_k-{\mathbf{V}}_j\oplus {\mathbf{m}}_k^{\prime }\Vert }^2=\sum _{i=1}^m{\Vert {\mathbf{Z}}_k({\mathbf{V}}_i,\hspace{0.17em}{\mathbf{V}}_j)\Vert }^2$$

(8)

Here, the m_k, m′_k is the pair of closest points in their respective local coordinate system, and N defines a small neighborhood of points taken in the order of hundreds of milliseconds that is assumed to have negligible pose error. Since we will require a linear approximation of this function and ⊕ is decidedly non-linear, we apply a Taylor expansion of Z_k:

$${\mathbf{Z}}_k({\mathbf{V}}_i,\hspace{0.17em}{\mathbf{V}}_j)\approx {\mathbf{Z}}_k({\overline{\mathbf{V}}}_i,\hspace{0.17em}{\overline{\mathbf{V}}}_j)-({\nabla }_{{\mathbf{V}}_i}{Z}_k({\overline{\mathbf{V}}}_i,\hspace{0.17em}{\overline{\mathbf{V}}}_j)\Delta {\mathbf{V}}_i-{\nabla }_{{\mathbf{V}}_i}\hspace{0.17em}{\mathbf{Z}}_k({\overline{\mathbf{V}}}_i,\hspace{0.17em}{\overline{\mathbf{V}}}_j)\Delta {\mathbf{V}}_j)$$

Here, ∇_Vi refers to the derivative with respect to V_i. We create a matrix decomposition ∇_Vi Z_k = M_k H, such that the matrix, M_k, is independent of V_i. Then, Z_k(V_i, V_j) is approximated as:

$${\mathbf{Z}}_k({\mathbf{V}}_i,\hspace{0.17em}{\mathbf{V}}_j)\approx {\mathbf{Z}}_k({\overline{\mathbf{V}}}_i,\hspace{0.17em}{\overline{\mathbf{V}}}_j)-{\mathbf{M}}_k({\mathbf{H}}_i\Delta {\mathbf{V}}_i-{\mathbf{H}}_j\Delta {\mathbf{V}}_j)$$

The minimum of this new linearized formulation of the error metric and its covariance is given by:

$$\begin{array}{l}{\overline{\mathbf{E}}}_{i,j}={({\mathbf{M}}^T\hspace{0.17em}\mathbf{M})}^{-1}{\mathbf{M}}^T\hspace{0.17em}\mathbf{Z}\\ {\mathbf{C}}_{i,j}=s^2({\mathbf{M}}^T\hspace{0.17em}\mathbf{M})\end{array}$$

Here, Z is the concatenated vector consisting of all Z_k = V̄_i ⊕ m_k – V̄_j ⊕ m^′_k. s is the unbiased estimate of the covariance of the independent and identically distributed errors of the error metric and is computed as:

$$s^2={(\mathbf{Z}-\mathbf{M}{\overline{\mathbf{V}}}_{i,j})}^T(\mathbf{Z}-\mathbf{M}{\overline{\mathbf{V}}}_{i,j})/(2M-3)$$

(9)

The matrix decomposition, M_k H, depends on the parametrization of the pose compounding operation ⊕. Since we employ the Euler parametrization as given in Equation (6), the matrices, M_k and H_i, are given by:

$$M_k=\left(\begin{array}{cccccc}1& 0& 0& 0& -m_{y,k}& -m_{z,k}\\ 0& 1& 0& m_{z,k}& m_{x,k}& 0\\ 0& 0& 1& -m_{y,k}& 0& m_{x,k}\end{array}\right)$$

(10)

$${\mathbf{H}}_i=\left(\begin{array}{cccccc}1& 0& 0& 0& {\overline{t}}_{z,i}\hspace{0.17em}\text{cos}({\overline{\theta }}_{x,i})+{\overline{t}}_{y,i}\hspace{0.17em}\text{sin}({\overline{\theta }}_{x,i})& {\overline{t}}_{y,i}\hspace{0.17em}\text{cos}({\overline{\theta }}_{x,i})\hspace{0.17em}\text{cos}({\overline{\theta }}_{y,i})-{\overline{t}}_{z,i}\hspace{0.17em}\text{cos}({\overline{\theta }}_{y,i})\hspace{0.17em}\text{sin}({\overline{\theta }}_{x,i})\\ 0& 1& 0& -{\overline{t}}_{z,i}& -{\overline{t}}_{x,i}\hspace{0.17em}\text{sin}({\overline{\theta }}_{x,i})& -{\overline{t}}_{x,i}\hspace{0.17em}\text{cos}({\overline{\theta }}_{x,i})\hspace{0.17em}\text{cos}({\overline{\theta }}_{y,i})-{\overline{t}}_{z,i}(\hspace{0.17em}\text{sin}{\overline{\theta }}_{y,i})\\ 0& 0& 1& {\overline{t}}_{y,i}& -{\overline{t}}_{x,i}\hspace{0.17em}\text{cos}({\overline{\theta }}_{x,i})& {\overline{t}}_{x,i}\hspace{0.17em}\text{cos}({\overline{\theta }}_{y,i})\hspace{0.17em}\text{sin}({\overline{\theta }}_{x,i})+{\overline{t}}_{y,i}\hspace{0.17em}\text{sin}({\overline{\theta }}_{y,i})\\ 0& 0& 0& 1& 0& \text{sin}({\overline{\theta }}_{y,i})\\ 0& 0& 0& 0& \text{sin}({\overline{\theta }}_{x,i})& \text{cos}({\overline{\theta }}_{x,i})\hspace{0.17em}\text{cos}({\overline{\theta }}_{y,i})\\ 0& 0& 0& 0& \text{cos}({\overline{\theta }}_{x,i})& -\text{cos}({\overline{\theta }}_{y,i})\hspace{0.17em}\text{sin}({\overline{\theta }}_{x,i})\end{array}\right)$$

(11)

4.2. Pose Optimization

We maximize the likelihood of all pose estimates, V_i,j, and their respective covariances, C_i,j, via the Mahalanobis distance:

$$W=\sum _i\sum _j({\overline{\mathbf{V}}}_{i,j}-({{\mathbf{V}}^{\prime }}_i-{{\mathbf{V}}^{\prime }}_j)){\mathbf{C}}_{i,j}^{-1}({\overline{\mathbf{V}}}_{i,j}-({{\mathbf{V}}^{\prime }}_i-{{\mathbf{V}}^{\prime }}_j))$$

or, with the incidence matrix, H, in matrix notation:

$$W={(\overline{\mathbf{V}}-\mathbf{HV})}^T\hspace{0.17em}{\mathbf{C}}^{-1}(\overline{\mathbf{V}}-\mathbf{HV})$$

(12)

This linear formulation of the problem has been made possible by the linearization of the pose difference equation in the previous section. The minimization of W is accomplished via solving the following linear equation system:

$$({\mathbf{H}}^T\hspace{0.17em}{\mathbf{C}}^{-1}\mathbf{H})\mathbf{V}={\mathbf{H}}^T\hspace{0.17em}{\mathbf{C}}^{-1}\overline{\mathbf{V}}$$

(13)

Computing the optimized trajectory is then reduced to inverting a positive definite matrix [29]. Since a significant portion of correspondences, i, j, are not considered, due to the lack of point correspondences, the matrix is sparse, so that we make use of sparse Cholesky decomposition [30]. A consequence of this formulation is that an arbitrary pose, V_i, must remain fixed, otherwise Equation (13) would be underdetermined. This fixed pose also incidentally defines the global coordinate system. Thus, for the algorithm, the optimization problem is decoupled from determination of global frame coordinates of the point cloud.

5. Experiments

To evaluate the newly proposed system for improving the scan quality, we have acquired a variety of data using multiple mobile laser scanners. An overview of the datasets is given in

Table 1. Overview of datasets used in the experiments. OBD, on-board diagnosis; IMU, inertial measurement unit; GNSS, Global Navigation Satellite System.

Environment	Sensors	Data	Experiment
Basement 1	Irma3D—Riegl VZ-400 laser scanner, Xsens MTi IMU, odometry of Volksbot RT3	4 short runs; 0.1 m/s; 3 to 7 million points each	evaluation of calibration procedure
Basement 2	Irma3D—Riegl VZ-400 laser scanner, Xsens MTi IMU, odometry of Volksbot RT3	1 short run; 0.2 m/2; 3 to 7 million points each	quantitative evaluation of semi-rigid optimization
Campus	VW Touran—Riegl VZ-400 laser scanner, Xsens MTi IMU, OBD device	500 m; 19 km/h; 12 million points	qualitative evaluation of semi-rigid optimization with very bad input data
Salzburg	Riegl—Riegl VZ-400, GPS RTK		qualitative evaluation of semi-rigid optimization algorithm
Dortmund 1	Optech Lynx Mobile Mapper; twin 2D laser scanners; Applanix POS L with integrated odometry, IMU and GNSS	circle around a park; unbroken trajectory; 150 m	qualitative evaluation of semi-rigid optimization
Dortmund 2	Optech Lynx Mobile Mapper; twin 2D laser scanners; Applanix POS L with integrated odometry, IMU and GNSS	two runs along the same road; 400 m	extensive quantitative analysis of semi-rigid optimization
Vienna	Riegl—Riegl VMX-450; Applanix POS L with integrated odometry, IMU and GNSS	two runs around a street corner; 300 m	comparison of calibration to state-of-the-art method; qualitative evaluation of semi-rigid optimization

. The first scanner is the autonomous robot platform Irma3D, as seen in Figure 2. The main sensor of Irma3D is a Riegl VZ-400 laser scanner. Irma3D is built with a Volksbot RT-3 chassis and uses the Xsens MTi IMU, as well as odometry to estimate its position. Several datasets using this robot were acquired in an empty basement room (see Figure 2). All datasets were acquired in continuous mode, i.e., the laser scanner rotates around its vertical axis while the robot moves simultaneously. The robot moved in “serpentine” trajectories, i.e., taking left and right turns, as well as segments where the heading remains unchanged. In this paper, we present five datasets recorded by Irma3D with three to seven million points each. Four of those datasets were acquired to demonstrate the calibration procedure as presented in this paper. For these datasets, the robot drove in a a relatively slow manner (approximately 0.1 m/s), to ensure minimal wheel slippage. The fifth dataset is used in the evaluation of the optimization algorithm. In this case, the robot was deliberately moving faster (approximately 0.2 m/s) to produce more errors in the pose estimation. For all five datasets, the scanner needed 6 s per rotation. 2D scans with a resolution of 0.1° were acquired at a rate of 120 Hz, which equates to a vertical resolution of 0.5° for the 3D point cloud.

The second mobile scanner uses the same sensors and scanning principle. However, as seen in Figure 3, the platform used was a Volkswagen (VW) Touran automobile. The vehicle is not equipped with GNSS or accurate odometry. The current forward velocity is acquired with an accuracy of only 1 km/h via an on-board diagnostics (OBD) device that is directly connected to the entertainment Controller Area Network (CAN bus)of the car. The vehicle moved alongside the campus buildings, made a U-turn, moved straight back on the other side of the road and made another U-turn to return to the origin. The car drove a distance of about 500 m and acquired a total of 12 million points with an average velocity of about 19 km/h. The scanner rotation and resolution settings were the same as in the Irma3D datasets. Due to the large inaccuracies in the initial dataset, a quantitative ground truth comparison is not applicable. Nonetheless, we acquired a ground truth dataset of the campus using a terrestrial laser scanner at twelve positions to provide a point of reference.

We used data collected with an experimental platform constructed by RIEGLLaser Measurement Systems GmbH that is similar to the previous mobile scanners. This platform also employs the Riegl VZ-400 laser scanner in continuous mode. However, unlike the previous platforms, the vehicle has no odometry at all. Since there is also no IMU installed, only a GPS Real Time Kinematic (RTK) system was used for the estimation of the vehicle trajectory. This platform acquired a dataset in Salzburg that we optimized using the semi-rigid matching algorithm.

In addition to these datasets, we evaluated the optimization algorithm on two datasets acquired by TopScan GmbH using the Optech Lynx Mobile Mapper (see Figure 4). The Lynx mobile mapper employs twin 2D laser scanners and pose estimation by an Applanix POS L with integrated odometry, IMU and GNSS. The mobile mapping system traversed two streets in the city of Dortmund. In the first dataset, the vehicle encircled a small park in an unbroken trajectory with a length of about 150 m. For the second, larger dataset, the vehicle traversed the same road twice. Both trajectories were stitched together on their eastern end to create a larger dataset with a trajectory of a length of 400 m.

Lastly, we also show the suitability of the calibration algorithm, as well as the optimization algorithm to the Riegl VMX-450. This system also uses an Applanix system for pose estimates with integrated odometry, IMU and GNSS. The data used for the experiments is a subset of a much larger dataset and consists of a single street corner that was traversed twice by the mobile laser scanner with a trajectory approximately 300-m in length. However, these two strips are sufficient to perform calibration of the mobile laser scanner using our algorithm.

5.1. Calibration Results

The calibration algorithm is evaluated with an indoor dataset, where a quantitative analysis is carried out, as well as on an outdoor dataset, where it is compared to a state-of-the-art calibration.

5.1.1. Calibrating the Mobile Robot Irma3D

We test the calibration algorithm on 4 indoor datasets acquired by Irma3D. A manual estimation of the calibration parameters had been performed before the experiments. We applied the automatic calibration technique to each dataset. A top view of two of the point clouds obtained with the automatically determined calibration parameters are shown in Figure 5. To evaluate the quality of the resulting point clouds, we compare them with a highly accurate model of the room as acquired by the Riegl VZ-400 laser scanner. The accuracy of the scanner is 5 mm [31] and, thus, the model should have a similar accuracy. Although the same scanner is used on the mobile robot, it is used in continuous mode, i.e., the laser scanner rotates around its vertical axis while the robot moves simultaneously. We compare the acquired point clouds with the model in the following fashion. After calibration, the entire mobile point cloud is matched to the model using ICP [20] from the 3D Toolkit (3DTK [32]). We compute point to plane distances on each of the 4 walls, as well as the ceiling and floor of the room. The deviations for single scans from the point cloud are plotted as color coded images, i.e., green for absolute errors less than 1 cm, yellow to red for large positive errors and cyan to blue for large negative errors. The full color scale is given in Figure 6. White areas indicate that no point was measured at the corresponding location.

The results of the direct comparison between representative scans from each of the four runs after automatic and manual calibration and the model of the room are shown in Figure 6. On the whole, the quality of the scans improves with the automatically determined parameters when compared to the manual estimation. Absolute errors are generally within 1–2 cm. Occasionally, the deviations exceed that boundary. Very rarely they are above 3–4 cm. These large-scale errors are explained by odometry related uncertainties in the estimation of the robot pose. Slipping wheels and shaking of the robot base are not detectable by any of the on-board sensors. Thus, errors in the pose estimation carry over to erroneous point clouds. The combined mean error for all datasets was improved by the automatic calibration from 10.3 mm to 8.3 mm.

5.1.2. Riegl VMX-450

We also demonstrate the proposed calibration algorithm on an outdoor dataset that was acquired using the Riegl VMX-450 in the city of Vienna. Figure 7 (bottom right) shows the point cloud as it was calibrated by RIEGLLaser Measurement Systems using state-of-the-art methods [6]. For this experiment, we consider these calibration parameters as the ground truth. We applied the calibration algorithm a total of 100 times to this dataset using ground truth parameters with additional random noise as the starting estimates. Noise was generated with a uniform distribution between −10° and 10° for each of the Euler angles that determine the orientation of the laser scanner. Other parameters, like the wheel circumference of the vehicle, could not be calibrated, as the necessary raw measurements were unfortunately not available for this dataset. The resulting calibration parameters are summarized in Figure 8 and

Table 2. Result of the outdoor boresight calibration experiment in degrees.

	Riegl in-House Calibration	Calibration Average	Calibration Median
θx	0.03286	0.01898	0.00878
θy	−0.14444	−0.13449	−0.09825
θz	0.12100	0.07695	0.08544

. Clearly, the proposed algorithm is capable of finding parameters very similar to those of the state-of-the-art algorithm with a high degree of certainty. Figures 7 and 9 show a representative starting estimate, the corresponding calibrated result, as well as the original ground truth dataset.

5.2. Optimization Results

We evaluate the algorithm using both indoor and outdoor datasets, where ground truth is available and a quantitative analysis of the produced point clouds is possible. We also provide other outdoor datasets, where no quantitative analysis is possible. Nonetheless, qualitative analysis on a variety of datasets shows the suitability of the optimization algorithm.

5.2.1. Trajectory Optimization of Spinning Scanners on Mobile Platforms

For a direct evaluation of the algorithm with comparison to ground truth data, we acquired another dataset (Basement 2) using the robot Irma3D in the empty basement room (see Figure 2). Ground truth data of the room was obtained using static terrestrial scanning with a Riegl VZ-400, with an estimated model accuracy of about 5 mm. To evaluate the quality of the resulting point cloud, we compare it to the high accuracy ground truth model of the room in the same fashion as in Section 5.1.1.

A qualitative comparison of the results of the algorithm is presented in Figure 10. The point clouds obtained with Irma3D along the full trajectory of dataset Basement 2 in the enclosed room are shown. The figure presents the initial point cloud (registered using robot IMU and odometry only), the rigid registration obtained with ICP, the semi-rigid registration computed by the novel algorithm and the ground truth data. The results of the direct comparison between the five individual point clouds, i.e., single rotations, from Basement 2 before and after automatic semi-rigid registration and the model of the room are shown in Figure 11. Again, the resulting images are color coded, i.e., green for absolute errors less than 1 cm, yellow to red for larger errors away from the observer and cyan to blue for larger errors towards the observer. White areas indicate that no point was measured at the corresponding location.

The proposed semi-rigid registration produces point clouds that more closely resemble the ground truth model. Although rigid registration procedures improve map quality to a degree, it is clear that the error within one scanner rotation is too large to allow for a good rigid match. The reduction of the inner error is even more obvious in the quantitative evaluation, as seen in Figure 11. The left side shows the point-to-model error of each of the individual scans before semi-rigid registration. The right column displays the error after the dataset has been corrected. Apart from a small region in the top scan, all errors have been reduced significantly by the semi-rigid registration.

In addition to the robots, we have data from two more mobile laser scanners that employ the VZ-400 in a continuously spinning mode. One is the automobile dataset, Campus, acquired on the campus of the Jacobs University Bremen. The initial and the semi-rigidly corrected point cloud, as well as the ground truth model of the campus, are shown in Figure 12. Due to the limited sensor capabilities of the mobile scanner, the initial trajectory is so erroneous that any attempts to improve upon the map quality by rigid registration methods failed. The semi-rigid registration shows a drastic improvement of the input data. The resulting point cloud exhibits no sign of error accumulation and reflects the ground truth quite well. It is remarkable that the algorithm also produces a straight trajectory when the only input is the first curved half of the trajectory. The scanner sweeps the same environment multiple times in 6 s intervals, such that local improvements in the trajectory are possible. The algorithm is capable of computing a valid trajectory with an input that is erroneous all the time.

Finally, the VZ-400 is also used on an experimental platform by the RIEGLLaser Measurement Systems GmbH equipped only with a GPS unit for trajectory estimation. An overview of the dataset is given in Figure 13. It is clear that the semi-rigid optimization algorithm improves on the quality of the point cloud. Figure 14 gives a more detailed view of the scene. Although the final registration is not perfect, the errors in the initial point cloud are reduced. The semi-rigid registration does not further improve this dataset, due to being trapped in a local minimum. This is an inherent limitation of the proposed algorithm, similar to the behavior of classical rigid registration using the ICP algorithm.

5.2.2. Lynx Mobile Mapper Data

The optimization algorithm was applied to the Dortmund datasets acquired by TopScan GmbH using the Optech Lynx Mobile Mapper. These datasets differ significantly from previous ones in that the Lynx Mobile Mapper employs two 2D laser scanners and the 3D point cloud is created purely by the movement of the vehicle. Datasets like this are produced by current state-of-the-art mobile scanning systems. Standard rigid optimizations generally do not apply, because there is only a small overlap between subsequent “sub-scans”. For these datasets, no information about the measurement certainty of the vehicle pose was given. The covariances, C_i,i+1, were set to a manually estimated diagonal matrix that reflect the estimated certainty of the pose estimation in general. However, the semi-rigid optimization proposed in this paper can deal with these challenging datasets, as well.

We first evaluate the algorithm on the smaller loop around the small park (Dortmund 1). Here, no ground truth data is available. However, in a direct comparison between the acquired data and the optimized data, one can see many improvements in scan quality; see Figure 15 for an overview. The measured scene was static for the observational period, i.e., the parked cars did not move. The most noticeable improvements are in features like tree trunks and advertising columns, which should have a circular shape in the orthogonal projection in Figure 15. Comparing building facades before and after optimization also reveals an improvement in scan quality. The video supplement to this paper is an animated fly-through of the dataset. It can be seen in high quality at [33].

For the Dortmund 2 dataset, the mobile scanner traversed the same scene twice in a straight trajectory. This minimizes errors due to rotational uncertainties and is a more realistic use of the mobile scanning system. An overview and closeup of the dataset before and after semi-rigid optimization is given in Figure 16. For this dataset, a quantitative analysis is available. Planes were extracted for the facades of the buildings. This was done for both of the two 2D scanners, as well as both of the passes. Thus, every facade is represented up to four times in the dataset. Comparing two plane parameters of the same facade yields an error vector with a direction and length and gives a measure of relative accuracy. With four planes, there can be a total of six non-redundant comparisons. From these, we exclude the two comparisons between planes extracted from different sensors, but in the same pass. Error vectors between these reflect only how good the system is calibrated. Compared to the other evaluations, the errors are very small and are not improved by the semi-rigid optimization. This analysis has been conducted by TopScan GmbH. The result of this analysis for the initial dataset is given in Figure 17. The average distance error between the facades is 10.6 cm. Averages and deviations for each of the four comparisons is given in Figure 18. The same analysis was carried out on the optimized dataset. The results are given in Figure 19. From a direct comparison of the results, it is clear that the optimization succeeded in improving the quality of the point cloud. Averages and variances are also given in Figure 18. The average error is reduced to only 2.3 cm. The standard deviation is reduced from 11.5 cm to 2.6 cm. Even the maximum errors are reduced significantly. Similar to the indoor datasets, we have also plotted the deviation between the first and second pass of the facades of the nearby buildings. The deviations both before and after optimization are color coded and displayed in Figure 20.

In addition to the evaluation of the inner errors of the point clouds, we also evaluated their external accuracy. This was done by way of five control points that were distributed in the scene. Each of the control points were located within a global reference with an absolute accuracy of approximately 1 cm by a professional surveyor on site. The control points were also manually located in the point clouds of both 2D scanners for both passes. This is with the exception of a single control point that was not seen by one scanner in one of the passes. Thus, a total of 19 distance errors for the original dataset, as well as the optimized dataset, were determined. The result of this evaluation is given in

Table 3. Differences between the position of the control points on site and the position in the point cloud, both for the original Lynx Mobile Mapper data (Left) and the optimized point cloud (Right).

	Original	Optimized
L1–L2	0.132	0.178
	0.121	0.113
	0.118	0.171
	0.238	0.267
	0.245	0.171
L1–R2	0.121	0.171
	0.130	0.115
	0.105	0.156
	0.241	0.267
	0.282	0.202
R1–L2	0.099	0.124
	0.139	0.123
	0.151	0.169
	0.214	0.254
	0.156	0.158
R1–R2	0.099	0.129
	0.147	0.130
	0.226	0.242
	0.157	0.167
average standard deviation	0.164	0.174
	0.055	0.049

. The difference between the original and the optimized dataset has a minimum of less than 1 cm, while the variances are 5.6 cm or 4.9 cm, respectively. Given the small number of independent sample points, we can conclude that the global error did not change to a statistically significant degree. This behavior is encouraging, since the optimization algorithm did not currently take into account the certainty with which the pose of the vehicle was measured. Furthermore, the coordinate system was fixed with respect to an arbitrary pose estimate, in our case the beginning of the trajectory of the first pass. We expect that the algorithm will also improve on the external accuracy of a given dataset when GNSS measurements are used to define the global coordinate system and GNSS measurement certainty is taken into account.

5.2.3. Riegl VMX-450

We also applied the semi-rigid optimization to data that was acquired by the Riegl VMX-450 laser scanner system. This Vienna dataset is depicted in Figure 21. Unfortunately, no quantitative evaluation of this dataset is available. However, as in previous experiments, we can visually inspect the results of the semi-rigid registration. This particular dataset is initially of relative high quality, but some errors are still present. Applying the optimization algorithm reduces the occurrence of surfaces appearing multiple times and leads to a more appealing point cloud.

6. Conclusions

The fully automatic calibration algorithm for mobile laser scanners as presented in this paper has proven to increase point cloud quality. It is capable of finding more precise calibration parameters than manual estimation. The algorithm is robust against non-systematic errors that are necessarily present in any mobile laser scanning system, such as laser scanner noise and erroneous pose estimations due to slipping wheels. A minor issue of the algorithm is the dependency on the environment. In our experiments, the robot primarily moved on a level surface. This means some calibration parameters, i.e., the upwards translation of the laser scanner in relation to the vehicle frame, are less constrained by the quality measure as other parameters.

The formulation and implementation of the algorithm is generalized to allow for a wide range of sensors to be calibrated. In future work, we will study how well the algorithm deals with multiple sensors and additional sensors of different modality.

The proposed semi-rigid registration algorithm has shown that it is capable of processing and significantly improving upon a variety of datasets. It exceeds current state-of-the-art rigid registration techniques, not only when comparing the final maps produced, but also when comparing the inner deformation of subsections of the point cloud with ground truth. In future work, we will incorporate information about pose measurement certainty, especially GNSS error estimates, to improve the overall global accuracy. Since the semi-rigid registration may converge only to a local minimum, the conditions for the convergence must be inspected in more detail in the future. Furthermore, we plan to compare the precise maximum likelihood 3D point clouds with some reference independent points using a global coordinate system and aim at developing algorithms for improving the global accuracy of the point clouds.