Automatic Multi-Sensor Extrinsic Calibration for Mobile Robots

TL;DR

This paper introduces an automatic method for calibrating multiple heterogeneous sensors on mobile robots, estimating full 6DoF parameters efficiently and accurately, even with monocular cameras and in various environments.

Contribution

It presents a novel motion-based calibration approach that initializes parameters in closed form and explicitly handles scale ambiguity, enabling multi-sensor fusion on mobile robots.

Findings

Method achieves accurate calibration in simulation and real-world tests.

Outperforms existing motion-based calibration methods.

Handles monocular camera scale ambiguity explicitly.

Abstract

In order to fuse measurements from multiple sensors mounted on a mobile robot, it is needed to express them in a common reference system through their relative spatial transformations. In this paper, we present a method to estimate the full 6DoF extrinsic calibration parameters of multiple heterogeneous sensors (Lidars, Depth and RGB cameras) suitable for automatic execution on a mobile robot. Our method computes the 2D calibration parameters (x, y, yaw) through a motion-based approach, while for the remaining 3 parameters (z, pitch, roll) it requires the observation of the ground plane for a short period of time. What set this proposal apart from others is that: i) all calibration parameters are initialized in closed form, and ii) the scale ambiguity inherent to motion estimation from a monocular camera is explicitly handled, enabling the combination of these sensors and metric ones (Lidars, stereo rigs, etc.) within the same optimization framework. %Additionally, outlier observations arising from local sensor drift are automatically detected and removed from the calibration process. We provide a formal definition of the problem, as well as of the contributed method, for which a C++ implementation has been made publicly available. The suitability of the method has been assessed in simulation an with real data from indoor and outdoor scenarios. Finally, improvements over state-of-the-art motion-based calibration proposals are shown through experimental evaluation.

Tables

Table 1. TABLE I: Estimated calibration parameters w.r.t. the odometer

Noise		$x_x$	$x_y$	$x_z$	$x_{\theta }$	$x_{\beta }$	$x_{\alpha }$	$x_s$
Ground truth		$50\mathrm{cm}$	$10\mathrm{cm}$	$100\mathrm{cm}$	$-90°$	$4.77°$	$-135°$	$2$
$\lambda =0$	Avg.	$50\mathrm{cm}$	$10\mathrm{cm}$	$100\mathrm{cm}$	$-90°$	$4.77°$	$-135°$	2
	[8]	$44.6\mathrm{cm}$	$10\mathrm{cm}$	–	$-95.2°$	$-12.2°$	$-150.9°$	–
	Avg.	$49.2\mathrm{cm}$	$9.8\mathrm{cm}$	$99.5\mathrm{cm}$	$-89.8°$	$4.77°$	$-135°$	$1.99$
$\lambda =1$	Abs.	$1.0\mathrm{cm}$	$0.2\mathrm{cm}$	$0.5\mathrm{cm}$	$0.5°$	$0.0°$	$0.01°$	$0.01$
	Rel.	$2.15\%$	$2.45\%$	$0.6\%$	$0.6\%$	$0.0\%$	$0.01\%$	$0.6\%$
	Avg.	$46.8\mathrm{cm}$	$9.4\mathrm{cm}$	$98.4\mathrm{cm}$	$-90.27°$	$4.77°$	$-134.9°$	$1.97$
$\lambda =2$	Abs.	$3.4\mathrm{cm}$	$0.7\mathrm{cm}$	$1.6\mathrm{cm}$	$0.7°$	$0.0°$	$0.04°$	$0.03$
	Rel.	$6.9\%$	$7.1\%$	$1.6\%$	$0.8\%$	$0.0\%$	$0.03\%$	$1.5\%$

Table 2. TABLE II: Estimated calibration parameters w.r.t. the vehicle (translations are in  m meter \mathrm{m} and rotations in degrees)

Calibration		Camera
		$x_x$	$x_y$	$x_z$	$x_{\theta }$	$x_{\beta }$	$x_{\alpha }$
Ground truth		$2.21$	$0.43$	$2.25$	$-88.4$	$-3.0$	$-87.2$
	parkin_0L	$2.19$	$0.43$	$2.37$	$-88.4$	$-3.7$	$-87.8$
	parkin_2L	$2.20$	$0.43$	$2.37$	$-88.4$	$-3.7$	$-87.8$
Ours	parkin_6L	$2.19$	$0.43$	$2.37$	$-88.4$	$-3.7$	$-87.8$
	Avg.	$2.20$	$0.43$	$2.37$	$-88.4$	$-3.7$	$-87.8$
	Stdv.	$0.01$	$0.00$	–	$0.1$	–	–
	parkin_0L	$2.20$	$0.42$	–	$-87.5$	$-2.8$	$-87.3$
	parkin_2L	$2.21$	$0.39$	–	$-90.3$	$-3.0$	$-87.1$
Della Corte	parkin_6L	$2.22$	$0.37$	–	$-93.2$	$-0.7$	$-87.0$
et al. [8]	Avg.	$2.21$	$0.39$	–	$-89.5$	$-3.0$	$-87.1$
	Stdv.	$0.01$	$0.03$	–	$1.76$	$0.21$	$0.15$

Table 3. TABLE III: Estimated calibration parameters w.r.t. the odometer (translations are in  cm centimeter \mathrm{cm} and rotations in degrees)

Calibration		Laser scanner			Monocular Camera
		$x_x$	$x_y$	$x_{\theta }$	$x_x$	$x_y$	$x_z$	$x_{\theta }$	$x_{\beta }$	$x_{\alpha }$
Expected		$0.0$	$0.0$	$0.0$	$10.0$	$0.5$	$175$	$-95$	$2.5$	$-130$
	#1	$-0.4$	$1.2$	$-3.0$	$11.1$	$0.5$	$175.6$	$-93.7$	$3.85$	$-125.4$
	#2	$-0.3$	$1.7$	$-2.3$	$11.0$	$0.2$	$175.6$	$-93.6$	$3.85$	$-125.4$
Ours	#3	$0.0$	$1.1$	$-2.9$	$11.6$	$0.4$	$175.6$	$-93.9$	$3.85$	$-125.4$
	Avg.	$-0.2$	$1.3$	$-2.8$	$11.3$	$0.4$	$175.6$	$-93.8$	$3.85$	$-125.4$
	Stdv.	$0.2$	$0.3$	$0.4$	$0.34$	$0.15$	–	$0.13$	–	–
[8]	Avg.	$-2.8$	$1.4$	$-2.9$	$8.4$	$0.4$	–	$-101.6$	$-17.2$	$-137.3$
	Stdv.	$0.2$	$0.5$	$0.6$	$0.33$	$0.24$	–	$0.5$	$0.07$	$0.08$

Equations

$${}^i\mathbf{m}=x_s\left(R(x_{\theta }){}^j\mathbf{m}+\left[\begin{array}{c}{x}_x\\ x_y\end{array}\right]\right),$$

$$R(\theta )\triangleq \left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right].$$

$$\begin{array}{rl}\oplus :& \text{Sim(2)}\times \text{Sim(2)}\to \text{Sim(2)}\\ & \mathbf{a}\oplus \mathbf{b}=\left[\begin{array}{c}{b_s}^{-1}{a}_x+b_x\cos (a_{\theta })-b_y\sin (a_{\theta })\\ {b_s}^{-1}{a}_y+b_x\sin (a_{\theta })+b_y\cos (a_{\theta })\\ a_{\theta }+b_{\theta }\\ a_s{b}_s\end{array}\right]\end{array},$$

$$\begin{array}{rl}\ominus :& \text{Sim(2)}\to \text{Sim(2)}\\ & \ominus \mathbf{a}=\left[\begin{array}{c}-a_s{a}_x\cos (a_{\theta })-a_s{a}_y\sin (a_{\theta })\\ a_s{a}_x\sin (a_{\theta })-a_s{a}_y\cos (a_{\theta })\\ -a_{\theta }\\ {a_s}^{-1}\end{array}\right].\end{array}$$

$${\mathbf{p}}_i^k\triangleq \ominus {\mathbf{q}}_i^k\oplus {\mathbf{q}}_i^{k+1}.$$

$${\mathbf{q}}_j^k={\mathbf{q}}_i^k\oplus {}^i{\mathbf{x}}_j$$

$${\mathbf{p}}_j^k=\ominus \left({\mathbf{q}}_i^k\oplus {}^i{\mathbf{x}}_j\right)\oplus \left({\mathbf{q}}_i^{k+1}\oplus {}^i{\mathbf{x}}_j\right).$$

$${\mathbf{p}}_j^k=\ominus {}^i{\mathbf{x}}_j\oplus {\mathbf{p}}_i^k\oplus {}^i{\mathbf{x}}_j.$$

$${}^i{\mathbf{x}}_j^{\ast }=\underset{\mathbf{x}}{\mathrm{argmin}}{\mathcal{C}}_{ij}(\mathbf{x})$$

$${\mathcal{C}}_{ij}(\mathbf{x})\triangleq \frac{1}{2}\sum _k{\Vert {\mathbf{ϵ}}_{ij}^k(\mathbf{x})\Vert }_2^2$$

$${\mathbf{ϵ}}_{ij}^k(\mathbf{x})\triangleq {\mathbf{p}}_j^k-\ominus \mathbf{x}\oplus {\mathbf{p}}_i^k\oplus \mathbf{x},$$

$${\mathbf{R}}_i\triangleq R_y(x_{\beta })R_x(x_{\alpha }),{\mathbf{t}}_i\triangleq {\left[\begin{array}{ccc}0& 0& x_z\end{array}\right]}^{\top }$$

$${}^p{\mathbf{m}}_j={\mathbf{R}}_i{}^i{\mathbf{m}}_j+{\mathbf{t}}_i$$

$${}^i{\mathbf{x}}^{\ast }=\underset{\mathbf{x}}{\mathrm{argmin}}{\mathcal{E}}_i(\mathbf{x})$$

$${\mathcal{E}}_i(\mathbf{x})\triangleq \frac{1}{2}\sum _j{w}_j{\Vert {\mathbf{\eta }}_i^j(\mathbf{x})\Vert }_2^2,{\mathbf{\eta }}_i^j(\mathbf{x})\triangleq \mathbf{n}\cdot {}^p{\mathbf{p}}_j-D$$

$${\bar{\mathbf{R}}}_i^k\triangleq {\mathbf{R}}_i{\mathbf{R}}_i^k{\mathbf{R}}_i^{\top },{\bar{\mathbf{t}}}_i^k\triangleq {\mathbf{R}}_i{\mathbf{t}}_i^k,$$

$${\mathbf{ϵ}}_{ij}^k=\mathbf{x}\oplus {\mathbf{p}}_j^k-{\mathbf{p}}_i^k\oplus \mathbf{x}$$

$$\mathbf{\phi }={\left[\begin{array}{ccccc}{x_s}^{-1}& x_x& x_y& \cos x_{\theta }& \sin x_{\theta }\end{array}\right]}^{\top }$$

$${\mathbf{ϵ}}_k={\mathbf{Q}}_k\mathbf{\phi },$$

$${\mathbf{Q}}_k\triangleq \left[\begin{array}{ccccc}-{p_i^k}_x& 1-\cos {p_i^k}_{\theta }& \sin {p_i^k}_{\theta }& {p_j^k}_x& -{p_j^k}_y\\ -{p_i^k}_y& -\sin {p_i^k}_{\theta }& 1-\cos {p_i^k}_{\theta }& {p_j^k}_y& {p_j^k}_x\end{array}\right].$$

$${\mathcal{C}}_{ij}(\mathbf{x})=\frac{1}{2}{\mathbf{\phi }}^{\top }\mathbf{M}\mathbf{\phi }+C,$$

$$\mathbf{M}\triangleq \sum _k{{\mathbf{Q}}_k}^{\top }{\mathbf{Q}}_k.$$

$${\mathbf{\phi }}^{\ast }=\underset{\mathbf{\phi }}{\mathrm{argmin}}{\mathbf{\phi }}^{\top }\mathbf{M}\mathbf{\phi },\text{subject to}{\phi }_4^2+{\phi }_5^2=1.$$

$${\mathbf{\phi }}^{\top }\mathbf{W}\mathbf{\phi }=1,\text{with }\mathbf{W}\triangleq \left[\begin{array}{cc}{0}_{3\times 3}& 0_{3\times 2}\\ 0_{2\times 3}& {\mathbf{I}}_{2\times 2}\end{array}\right].$$

$$\mathcal{L}(\mathbf{\phi },\lambda )\triangleq {\mathbf{\phi }}^{\top }\mathbf{M}\mathbf{\phi }+\lambda \left({\mathbf{\phi }}^{\top }\mathbf{W}\mathbf{\phi }-1\right)$$

$$\frac{\partial \mathcal{L}}{\partial \mathbf{\phi }}=2{\mathbf{\phi }}^{\top }\left(\mathbf{M}+\lambda \mathbf{W}\right)=0^{\top }.$$

$$\det (\mathbf{M}+\lambda \mathbf{W})=0$$

$$\bm{\bar{\varphi}}^{*}=\frac{\operatorname*{sign}(\gamma_{1}^{*})}{\left\lVert\bigl{[}\begin{smallmatrix}\gamma_{4}^{*}&\gamma_{5}^{*}\end{smallmatrix}\bigr{]}^{\top}\right\rVert_{2}}\bm{\gamma}^{*}.$$

$${}^i{\mathbf{x}}_j^{\ast }={\left[\begin{array}{cccc}{\phi }_2^{\ast }& {\phi }_3^{\ast }& \text{atan2}({\phi }_5^{\ast },{\phi }_4^{\ast })& {{\phi }_1^{\ast }}^{-1}\end{array}\right]}^{\top }.$$

$$\mathbf{n}\triangleq {\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^{\top },D\triangleq 0$$

$${\mathbf{\eta }}_i^j={{}^p{m}_j}_z={\mathbf{r}}_z{}^i{\mathbf{m}}_j+x_z$$

$${\mathbf{r}}_z\triangleq \left[\begin{array}{ccc}-\sin x_{\beta }& \cos (x_{\beta })\sin (x_{\alpha })& \cos (x_{\beta })\cos (x_{\alpha })\end{array}\right]$$

Full text

I Introduction

Autonomous mobile robots require sensors to perceive the environment and estimate their state. Combining measurements from different sensors is a need to improve robustness and to compensate for individual sensor limitations. In order to express measurements into a common reference system, accurate relative transformations between the sensors are required (i.e. extrinsic calibration). For example, when doing Visual Odometry [1] it is a common practice to consider also information from an odometer [2, 3, 4], which requires knowing the spatial transformation between the camera and the wheel odometry system.

Two main approaches can be found in the literature for extrinsic calibration: i) exploiting a priori knowledge about the environment or ii) relying on per-sensor motion estimates. In the first case, a priori information is typically in the form of scene patterns [5, 6, 7] or known landmarks [4]. Pattern-based calibration methods are usually developed for specific pairs of sensors. For example: a camera and a 2D laser scanner [5], or monocular and depth cameras [6]. Using specific calibration techniques to cover all sensors pairs in a multi-sensor configuration makes the calibration process complex and almost impractical. Additionally, relying on known landmarks for calibration prevents the automatic execution of the calibration process since they have to be placed in the environment beforehand. However, automatic calibration methods are interesting for mobile robots: they demand re-calibration capabilities for long-term operation since the relative transformation between their sensing devices can change due to crashes, vibrations, human intervention, etc [8, 9].

In contrast to pattern-based calibration methods, motion-based ones do not need to modify the environment and they can be used with any sensing devices, provided that ego-motion can be inferred from them. Monocular cameras are the exception, being excluded due to the scale ambiguity problem, which has to be considered explicitly. Another practical limitation arises from the planar movement of mobile robots, which prevents the estimation of the full 6DoF extrinsic calibration parameters solely form incremental motions, as shown in [2].

In this paper, we contribute a method to extrinsically calibrate multiple heterogeneous sensors such as odometers, visual and depth cameras, 2D and 3D laser scanners, etc., mounted on a robot moving on a planar ground (see Fig. 1). Our approach computes the 6 extrinsic parameters in two main steps: the 2D calibration parameters (x, y, yaw) are obtained from per-sensor incremental motions, while the remaining 3 parameters (z, pitch, roll) are computed from the observation of the ground plane. The ground plane has to be observed for a short period of time, and this requisite is fulfilled by most robotic sensor configurations. By decoupling the extrinsic calibration problem into two simpler sub-problems we are able to solve all parameters in closed form, thus removing the necessity to provide initial estimates. Another important aspect of our method is that it formulates the calibration problem taking into account the scale ambiguity in monocular cameras, which is estimated in the motion-based step. Thus, our approach allows to combine these commonly-used devices with other types of sensors in the same calibration framework. A C++ open source implementation of the contributed method is publicly available for download at https://github.com/dzunigan/robot_autocalibration.

Our proposal has been validated in a simulation experiment and assessed with real data from two different datasets: one in outdoor scenarios containing information from a monocular camera and a GPS [10], and another one gathered indoors using one of our mobile robotic platforms, equipped with an odometer, a monocular camera and a 2D laser scanner [11]. We also report improvements in terms of accuracy and consistency of the estimates over the state-of-the-art motion-based approach from Della Corte et al. [8].

II Related Work

Most extrinsic calibration methods found in the literature consider specific sensors and exploit some a priori information about the observed scene. For example, the work in [5] uses scene corners (orthogonal trihedrons) to calibrate a 2D laser scanner and a camera. The information provided by a spherical object is used in [6] to calibrate a camera and a depth sensor. The overlap requirement between the sensors’ field of view is relaxed in [7] to calibrate multiple depth cameras from the common observation of planar structures. In contrast to these methods, a priori knowledge is not required by motion-based calibration approaches, and hence can be used automatically in unstructured environments.

The extrinsic calibration from individual sensor egomotion was first considered, to the best of our knowledge, by Brookshire and Teller in [12]. They present an observability analysis of the 2D extrinsic calibration problem and a solution based on the iterative minimization of a least-squares cost function. Similarly, in [13], Censi et al. consider the calibration of both, a generic sensor (for which egomotion can be estimated) with respect to the robot reference system, and the intrinsic parameters of a differential drive configuration. They also analyze the observability of the calibration problem and propose a closed form solution to a least-squares formulation. In [14], Brookshire and Teller extend their previous formulation by considering the calibration of the 3D transformation parameters. Schneider et al. present in [15] an online extrinsic calibration approach based on the Unscented Kalman Filter.

Unlike these previous works, the authors in [2] address the scale ambiguity of monocular cameras proposing a closed form least-squares solution to the odometer-camera calibration problem based on incremental motions. Their approach does not require any a priori information about the environment, and thus allows for an automatic execution. Zienkiewicz and Davison consider in [3] a similar calibration problem. They propose a method based on the parameterization of the homography. The solution is found by minimizing the photometric errors induced by the homography arising from planar motions observing the ground plane.

All aforementioned works are limited to the calibration of pairs of sensors. In contrast, the work in [16] considers the calibration of multiple cameras with an odometer. Their pipeline estimates the intrinsic camera parameters as well as the extrinsic parameters with respect to an odometer (initialized from [2]), but it is designed to work only with imaging sensors. The work in [17] tackles the calibration of multiple sensors explicitly. Their formulation is based on the Gauss-Helmert model, where all motions have to be expressed in the same scale, an thus it is not suitable for cameras. Taylor and Nieto propose in [18] an algorithm to extrinsically calibrate multiple sensors from 3D motions, including cameras, by solving multiple hand-eye calibration problems. Recently, Della Corte et al. presented in [8] an extrinsic calibration framework from 3D motions, where each sensor’s time delay are also estimated. However, as shown in [2], the full 6DoF calibration parameters cannot be estimated solely from planar motions performed by mobile robots.

In contrast to previous approaches, the method presented in this paper estimates the 6DoF calibration parameters of multiple, heterogeneous sensors (including monocular cameras) mounted on a mobile robot performing planar motions. For that, it relies on incremental sensor motions as well as the observation of the ground plane. The a priori information requirement introduced by the calibration step based on the ground plane detection is soft enough to allow automatic execution with different sensors (e.g. monocular cameras or 3D laser scanners) without requiring overlapping field of view configurations. Additionally, all extrinsic calibration parameters are estimated in closed form, and hence initial guess values are not required. As a consequence, our method is suitable for automatic execution by mobile robots.

III Problem Formulation

III-A Motion-based Calibration Coplanar Sensors

First, we consider the motion-based extrinsic calibration of two coplanar sensors under the assumptions that: i) 2D synchronized incremental poses are available for each of them, and ii) the translation component may be expressed in different scales. Our objective is to estimate the fixed 2D transformation between the two sensors that better explains the difference in the observed incremental motions.

More formally, our aim is to estimate the parameters of the 2D similarity transformation between the $i$-th and the $j$-th sensors ${}^{i}\bm{x}_{j}=\left(x_{x},x_{y},x_{\theta},x_{s}\right)\in\text{Sim(2)}$, where $x_{x},x_{y}\in{\rm I\!R}$ are the translation components, $x_{\theta}\in{\rm I\!R}$ the rotation angle, $x_{s}\in{\rm I\!R}^{+}$ the scaling factor (see Fig. 2), and Sim(2) the group of orientation-preserving similarity transformations111Similarity transformations are rigid transformations followed by a scaling [19].. The calibration transformation expresses the measurements taken from the $j$-th sensor ${}^{j}\bm{m}\in{\rm I\!R}^{2}$ into the $i$-th reference frame as [19]:

[TABLE]

where $R(\theta)$ is the 2D rotation matrix defined by the angle $\theta$:

[TABLE]

From (1), we can define the group operator of Sim(2) as:

[TABLE]

and its inverse as:

[TABLE]

We consider incremental rigid body motions as input for the calibration process, expressed as $\bm{p}^{k}_{i}\in\text{SE(2)}$ between sample times $k$ and $k+1$. Notice that a SE(2) transformation is a particular case of Sim(2) with a scaling factor set to the identity, thus the previously defined operators for Sim(2) also hold for SE(2) transformations.

The incremental motions can be derived from the sensor poses $\bm{q}^{k}_{i},\bm{q}^{k+1}_{i}\in\text{SE(2)}$ as222Throughout this paper, we set to $\ominus$ a higher precedence than $\oplus$ to improve the readability of expressions involving Sim(2) operations.:

[TABLE]

The poses corresponding to the $j$-th sensor can be expressed in terms of the $i$-th sensor poses and the extrinsic calibration between them:

[TABLE]

and, therefore, its incremental motion is given by:

[TABLE]

Rearranging terms in (7) and applying (5), the relative motions $\bm{p}^{k}_{i},\bm{p}^{k}_{j}\in\text{SE(2)}$ for the $i$-th and $j$-th sensors at a given time step $k$ and the extrinsic calibration parameters ${}^{i}\bm{x}_{j}\in\text{Sim(2)}$ are related by:

[TABLE]

This relationship allows us to establish a cost function for a least-squares formulation of the two-sensor extrinsic calibration problem:

[TABLE]

where the error function in (11) refers to the relationship between the observed incremental motions and the calibration parameters in (8).

III-B Sensor Coplanarity Relaxation

Now, consider the case of a sensor in a generic configuration, under the assumption that the plane in which motions are performed can be observed. Our objective is to set a common reference system to estimate the relative pose parameters in order to project estimated motions into this plane, and hence enforcing the coplanarity constraint.

Formally, we want to estimate the relative rigid body motion $(\bm{R}_{i},\bm{t}_{i})\in\text{SE(3)}$ between the $i$-th sensor and the ground plane. For convenience, we set the plane’s local reference system to be at the projection of the sensor’s origin into the plane and with its $z$-axis parallel to the plane’s normal pointing upwards (see Fig. 1). Note that the in-plane rotation (around the $z$-axis) can be arbitrarily set as it will be calibrated latter on (as the $x_{\theta}$ parameter).

Thus, ignoring the relative in-plane position and orientation, the remaining 3DoF parameters ${}^{i}\bm{x}=\left(x_{z},x_{\alpha},x_{\beta}\right)\in{\rm I\!R}^{3}$ for each sensor are the target of our estimation at this step. The parameter $x_{z}\in{\rm I\!R}^{+}$ represents the perpendicular distance and $x_{\alpha},x_{\beta}\in{\rm I\!R}$ two rotation angles. More specifically, we define the relative SE(3) transformation to be:

[TABLE]

where $R_{x}(\cdot)$ and $R_{y}(\cdot)$ represent parameterized 3D rotation matrices along the $x$ and $y$ axes, respectively.

Hence, the $j$-th 3D point ${}^{i}\bm{m}_{j}\in{\rm I\!R}^{3}$ observed by the $i$-th sensor can be expressed in the local ground coordinates by applying the relative rigid body transformation as:

[TABLE]

Based on this relationship, we can define a cost function for a weighted least-squares formulation of the coplanarity relaxation problem:

[TABLE]

where each $\bm{\eta}^{j}_{i}$ term represents the perpendicular distance of the $j$-th point ${}^{p}\bm{m}_{j}$ to the ground plane and $w_{j}\in{\rm I\!R}^{+}$ the weight. The ground plane is defined through the unit normal vector $\bm{n}\in{\rm I\!R}^{3}$ and the distance to the origin $D\in{\rm I\!R}^{+}$ (Hessian normal form). For convenience, we set the ground plane parameters to be $\bm{n}\triangleq(0,0,1)$ and $D\triangleq 0$ (see Section IV-B).

The observed 3D incremental motions $(\bm{R}^{k}_{i},\bm{t}^{k}_{i})\in\text{SE(3)}$ for the $i$-th sensor can be projected to the ground plane using the estimated rotation matrix $\bm{R}_{i}$ from (14) as:

[TABLE]

and then the SE(2) incremental motions can be easily recovered as the $x$-$y$ translation component of $\bm{\bar{t}}^{k}_{i}$, and by extracting the rotation angle about the $z$-axis from $\bm{\bar{R}}^{k}_{i}$.

These three parameters allow us to enforce the coplanarity constraint required by the motion-based calibration described in Section III-A and, therefore, to estimate the full 6DoF extrinsic parameters.

IV Method Description

In the proposed method, we estimate the 6DoF calibration parameters in two main steps. First, for each sensor, we estimate the parameters of the relative rigid body transformation with respect to the ground plane. Then, we estimate the extrinsic 2D parameters for each sensor with respect to the reference one from coplanar incremental motions. More specifically, the calibration pipeline can be summarized as follows:

Acquire sensor data while the robot is moving. 2. 2.

Estimate per-sensor egomotion. 3. 3.

Estimate z, pitch and roll relative to the ground plane. 4. 4.

Project estimated trajectories to the ground plane and resample synchronous incremental motions. 5. 5.

Perform a motion-based calibration of the remaining parameters: x, y, and yaw. 6. 6.

Refine the initial 2D parameters (x, y, yaw) in a joint optimization framework.

The rest of this section is structured as follows. In Section IV-A we show how to solve the calibration of two coplanar sensors from incremental motions. Next, in Section IV-B, we show how the parameters relative to the ground plane can be solved. Finally, in Section IV-C we discuss practical aspects, covering: the sampling of synchronous incremental motions, considerations for robust estimation, the final refinement, and the metric estimation of z in the monocular case.

IV-A Closed Form Solution to the Motion-based Pair-wise Calibration Problem

In order to solve (9), we follow a similar approach as the one described in [13]. The solution comes by first reducing the least-squares formulation to a quadratic system with a quadratic constraint. The constrained optimization problem is then uniquely solved in closed form by using the method of Lagrange multipliers.

To reduce the calibration problem to a quadratic system, we first rewrite the error terms in (11) as:

[TABLE]

Next, we parameterize the calibration angle $x_{\theta}$ by two independent variables: $\cos x_{\theta}$ and $\sin x_{\theta}$. Grouping the unknown parameters into the vector $\bm{\varphi}\in{\rm I\!R}^{5}$

[TABLE]

we can write the error terms (17) as:

[TABLE]

where the matrix $\bm{Q}_{k}$ contains the known coefficients:

[TABLE]

The cost function (10) can be written compactly as:

[TABLE]

where $C\in{\rm I\!R}$ is a constant term and $\bm{M}$ the symmetric matrix:

[TABLE]

Therefore, we can write the least squares problem (9) equivalently as a quadratic system with a quadratic constraint:

[TABLE]

The constraint in (23) corresponds to $\cos^{2}x_{\theta}+\sin^{2}x_{\theta}=1$ and can be written in matrix form as:

[TABLE]

Thus, the constrained optimization problem (23) can be solved considering the Lagrangian:

[TABLE]

and its necessary conditions for optimality:

[TABLE]

Since the scale factor $\varphi_{1}\in{\rm I\!R}^{+}$ has to be a positive real number, the $5\times 5$ matrix $\bm{M}+\lambda\bm{W}$ must be singular to satisfy (26). Thus, we solve for the value of $\lambda\in{\rm I\!R}$ that makes:

[TABLE]

and then we find the solution $\bm{\varphi}^{*}$ in the kernel of $\bm{M}+\lambda\bm{W}$.

For the matrix $\bm{W}$, the expression (27) is a second order polynomial in $\lambda$ and, therefore, can be solved in closed form. The two candidate values $\lambda^{*}$ are computed as the real roots of the second-order polynomial. The rank of the $5\times 5$ matrix $\bm{M}+\lambda^{*}\bm{W}$ is then at most 4 by construction. Provided that at least two linearly independent incremental motions are observed and that the two sensors follow nondegenerated trajectories, the matrix has rank exactly 4 and hence a one-dimensional kernel. For detailed information about the observability analysis we refer the reader to [12, 13].

To obtain the solution $\bm{\bar{\varphi}}^{*}$ associated to each $\lambda^{*}$, consider any nonzero vector $\bm{\gamma}^{*}$ in the kernel of $\bm{M}+\lambda^{*}\bm{W}$. Then, we impose the constraint (23) and the fact that $\bar{\varphi}_{1}^{*}$ must be positive to uniquely identify $\bm{\bar{\varphi}}^{*}$ from the nullspace:

[TABLE]

Then, we choose the correct solution $\bm{\varphi}^{*}$ to (23) from the two $\bm{\bar{\varphi}}^{*}$ as the one that achieves the lowest cost according to (21). Finally, we recover the optimal calibration parameters ${}^{i}\bm{x}^{*}_{j}\in\text{Sim(2)}$ between the considered pair of sensors from $\bm{\varphi}^{*}$ as:

[TABLE]

IV-B Closed Form Solution to the Coplanarity Relaxation Problem

In order to solve (14), we again reduce the least-squares formulation to a quadratic system with a quadratic constraint. Then, the constrained optimization problem is uniquely solved in closed form by using the method of Lagrange multipliers.

To reduce the relaxation problem to a quadratic system, we first start by simplifying the error function (15). Recalling from the problem formulation, the Hessian normal parameters of the ground plane are:

[TABLE]

and thus, the perpendicular distance $\bm{\eta}^{j}_{i}$ of the $j$-th point ${}^{p}\bm{m}_{j}$ to the ground plane reduces to its third component:

[TABLE]

where $\bm{r}_{z}$ represents the third row of the rotation matrix $\bm{R}\in\text{SO(3)}$ as defined in (12), for the $i$-th sensor.

Therefore, we parameterize the $x_{\alpha}$ and $x_{\beta}$ angles by three independent variables: $-\sin x_{\beta}$, $\cos(x_{\beta})\sin(x_{\alpha})$ and $\cos(x_{\beta})\cos(x_{\alpha})$. Grouping the unknown parameters into the vector $\bm{\varphi}\in{\rm I\!R}^{4}$

[TABLE]

we can write the error terms in matrix form

[TABLE]

and, then, the cost function in (15) yields:

[TABLE]

with the symmetric matrix $\bm{M}\triangleq\sum_{j}w_{j}Q_{j}^{\top}Q_{j}$.

Therefore, we can write the least squares problem (14) as the following quadratic system with a quadratic constraint:

[TABLE]

The constraint in (36) corresponds to the orthogonality property of SO(3) matrices and can be written as:

[TABLE]

For this particular matrix $\bm{W}$, the necessary conditions for optimality are characterized by a third order polynomial in $\lambda$ and therefore can be solved in closed form, with a maximum of three different candidate values $\lambda^{*}$. And as before, each solution $\bm{\bar{\varphi}}^{*}$ associated to a $\lambda^{*}$ can be uniquely recovered from any nonzero vector in the kernel of $\bm{M}+\lambda^{*}\bm{{W}}$ by imposing the orthogonality constraint (36) and the fact that the perpendicular distance $\bar{\varphi}^{*}_{1}$ has to be positive.

The optimal, minimal solution $\bm{\varphi}^{*}$ to (36) is chosen among all candidates $\bm{\bar{\varphi}}^{*}$ by evaluating the cost function (15). Finally, the calibration parameters for the $i$-th sensor ${}^{i}\bm{x}\in{\rm I\!R}^{3}$ easily follow from $\bm{\varphi}^{*}$.

IV-C Practical Considerations

First, the motion-based calibration procedure as described in Section IV-A requires synchronous incremental motions. However, heterogeneous sensors in general are asynchronous and have different sampling rates. We overcome this problem by setting a sensor to provide the time reference and then, for the other sensors, we resample synchronous incremental motions computed from linearly interpolated poses of the planar trajectories.

Secondly, motion-based calibration methods rely on trajectories estimated by other means. However, motion estimation algorithms are subject to large errors arising from local tracking failure or drift. In order to limit the impact that erroneous observations have on the calibration, a Random Sample Consensus (RANSAC) [20] framework can be used with our closed form solution. A RANSAC scheme requires a predefined threshold on the error terms in order to discard outliers. However, the error function as defined in (11) mixes the rotation error and the translation error, which have different magnitudes. What’s more, the translation error is defined in the space of the second sensor (denoted by $j$ in (11)), and has an arbitrary scale in the case of a monocular camera. Therefore, we use a slightly different error function for the RANSAC outlier detection step:

[TABLE]

which is the translation error expressed in the metric space of the $i$-th sensor. The error terms defined in this way depend on the translation as well as on the rotation parameters and allow us to set a predefined threshold intuitively, with a simple interpretation. Note that the error function defined in (11) is still used for the closed form solution, and (38) is used only for the inlier-outlier classification in a RANSAC framework.

Another important consideration for practical calibration of multiple sensors is the consistency of the transformations between all sensors. The motion-based calibration presented in Section IV-A only considers constraints between a sensor and the reference one. We can improve the consistency of the calibration by also considering constraints between the other sensors in a joint optimization framework.

Let the reference sensor have index 0. Assuming we consider $n$ sensors for calibration, we want to estimate the calibration parameters ${}^{0}\bm{x}_{1},\dots^{0}\bm{x}_{n}\in\text{Sim(2)}$. We define the joint calibration problem as:

[TABLE]

where $S$ represents the set containing all sensor pairs for which additional constraints are considered. The modified error function in (38) is used again to express the error terms of all sensors in the same metric space. Note that only pairs of sensors that have their first sensor providing metrically accurate motions may be considered for additional constraints. The expression $\ominus\,^{0}\bm{x}_{i}\oplus^{0}\!\bm{x}_{j}$ is just the relative transformation ${}^{i}\bm{x}_{j}$ expressed in terms of the calibration parameters ${}^{0}\bm{x}_{i},^{0}\bm{x}_{j}$. In our implementation, we used the Cauchy loss function $\rho$ to cope with unmodeled errors not detected during the RANSAC step. The joint calibration problem in (39) can be solved iteratively, starting from the closed form solution described in Section IV-A. For more information on how to solve this optimization problem, we refer the interested reader to [21].

Finally, in the monocular camera case, the metric value of $x_{z}$ can be recovered by applying the scale factor estimated in the motion-based calibration to the $x_{z}$ value observed during the coplanarity relaxation. However, the same reconstruction has to be used for both problems in order to guarantee that they share the same scale.

V Experimental Evaluation

We have conducted several experiments in order to demonstrate the suitability of our calibration method in practice. First, we validated our method with simulated data, where the ground truth parameters are known (Section V-A). Next, we assessed the calibration accuracy of our method in both outdoor (Section V-B) and indoor (Section V-C) environments with real data, and compared it with the state-of-the-art motion-based calibration approach proposed by Della Corte et al. [8].

V-A Validation with synthetic data

The goal of the simulation experiment is to characterize the accuracy of our method. Here, the calibration parameters are known and the sensors are perfectly synchronized. We considered a two-sensor setting for the validation: an odometer and a monocular camera. While the odometer provides metrically accurate motions, the camera yields motions in a different scale (imitating the scale ambiguity problem).

In the simulation, we commanded the robot to follow an eight-shaped path several times. The incremental motions are affected by unbiased, normally distributed noise in each axis independently. We considered different levels of noise by varying the standard deviation: $\lambda$0.1\text{,}\mathrm{cm} for translation and $\lambda$0.03\text{\,}\mathrm{rad} for rotation, given a noise level $\lambda\in{\rm I\!R}^{+}$. The odometer is affected by noise only in the $x$-$y$ axes for translation and in the $z$ axis for rotation, since it provides 2D motion estimates. On the other hand, the camera motion has noise in the three axes (for translation and rotation).

The simulated camera has a QVGA resolution and a diagonal FOV of $70.1\text{\,}\mathrm{\SIUnitSymbolDegree}$ with square pixels. For the coplanarity relaxation, the camera estimates a dense, scaled reconstruction of the ground plane in the same scale as the incremental motions. We added unbiased, normally distributed noise with $\lambda$1\text{,}\mathrm{cm}$$ standard deviation to the depth measurements (for a noise level $\lambda$).

We considered 10 independent runs for different levels of noise. The relative transformation between the sensors is fixed and shared by all runs. The estimated parameters, as well as the ground truth are listed in Table I. We also reported the average and the absolute and relative Root Mean Squared Error (RMSE). We included the noise-free case ($\lambda=0$) to show the correctness of our formulation and compared it with the method in [8]. The approach in [8] incorrectly detects a time delay and bias its solution. Additionally, it cannot handle higher levels of noise. On the other hand, our method retrieves the correct solution in the noise-free case and is more robust under noisy observations. From the table we can also conclude that the calibration parameters estimated from the coplanarity relaxation are more accurate than the motion-based ones. This is mainly due to a higher number of observations ($76\,800$ vs $74$). Finally, the calibration errors are concentrated in the 2D translation parameters, since they depend on the estimation of all the other parameters (apart from $x_{z}$).

V-B Outdoor evaluation

For the outdoor experimental evaluation, we chose the publicly available dataset in Blanco et al. [10], as it provides highly accurate trajectories for different rigidly mounted sensors (including cameras) and ground truth extrinsic calibrations. The dataset was recorded with an electric buggy while driving in both planar and nonplanar surfaces. The vehicle’s motion was estimated from a total of three Real Time Kinematics (RTK) GPS receivers and the ground truth extrinsic calibrations were initialized from manual measurements and then refined via nonlinear optimization techniques.

We selected the vehicle’s trajectory as the reference motion and used the trajectory of the left camera as provided by the dataset for extrinsic calibration. Note that even though the sensor being calibrated is a camera, we are using the trajectory provided by the dataset, which is metrically accurate. The vehicle’s 6DoF poses were recorded at $1\text{\,}\mathrm{Hz}$, while the images were captured with 1024$\times$768\text{\,}\mathrm{px} resolution at $7.5\text{\,}\mathrm{Hz}$.

As our method requires synchronous incremental motions, we synchronized the trajectory of the camera with the vehicle’s trajectory by linear interpolation. We chose such a trajectory as the synchronization time reference as it has the lower rate of the two and thus interpolated camera poses can be computed more accurately. Both the proposed method and the motion-based approach in [8] are given the same three input trajectories for the sake of fairness.

The calibration results of both methods for the three planar sequences are presented in Table II. For our method, the 3DoF calibration parameters of the relative transformation to the ground plane are estimated only once for all sequences (that is why no standard deviation values for these estimations are shown in the table). The 3D ground points are extracted from images by first running the Structure-from-Motion (SfM) pipeline [22] on 50 of them with enough texture on the ground and then detecting planar points on the reconstructed 3D scene through homographies. We used the metrically accurate trajectory of the camera to compute the scale factor in order to retrieve the metric value of $x_{z}$. That parameter is missing for [8], since it is unobservable solely from planar motions [2]. In general, both methods provide consistent parameters estimation. However, the method in [8] required a close initial guess of the translation component in order to provide reasonable results.

V-C Indoor evaluation

For the indoor experimental evaluation, we used the Giraff mobile robot [23, 24] equipped with an odometer (by integrating wheel encoders for a differential drive configuration), an Hokuyo UTW-30LX 2D laser scanner and an uEye UI-1240SE-M-GL monocular camera. The laser scanner is mounted in parallel to the ground, while the camera is pointing downwards with an incidence angle of about $40\text{\,}\mathrm{\SIUnitSymbolDegree}$, as depicted in Fig. 3 (left). The laser scanner and the camera are rigidly connected to the robot and, therefore, the relative transformation between them does not change during the experiments. Both the odometer and the camera have been intrinsically calibrated before attempting extrinsic calibration.

We considered five independent calibrations and then we analyzed the mean and deviation of the estimated parameters. To this end, we recorded three data sequences with software acquisition timestamps, while driving the robot through our lab following an eight-shaped path several times. The incremental motions of the laser scanner are estimated using the method in [25], which provides 2D pose estimates at about $6\text{\,}\mathrm{Hz}$. For the camera, we run ORB-SLAM [26] with 1280$\times$1024\text{\,}\mathrm{px} images at $5\text{\,}\mathrm{Hz}$.

This time we synchronized the incremental motions of the laser with respect to the keyframes as selected by the SLAM solution. We chose the camera as the synchronization time reference in order to interpolate between odometry and laser poses, which are available at higher and more consistent frame rates. Each data sequence contains $7090$ synchronized incremental motions.

The experimental results within indoor environments for the considered methods are presented in Table III, as well as the expected parameters for reference (from manual measurements). The robotic platform uses the laser for localization and thus contains rough calibration parameters for it. The motion estimation algorithm for the laser uses them to estimate the pose of the robot instead of the laser one. Thus, even though there is a translation of $22\text{\,}\mathrm{cm}$ in the $x$ axis, the expected calibration parameter is $x_{x}=$0\text{,}\mathrm{cm}$$. The $x_{z}$, $x_{\beta}$ and $x_{\alpha}$ parameters of the laser are omitted in Table III, since they are unobservable with the current setup [14]. Therefore, in the case of the laser scanner, we skipped the coplanarity relaxation for our method, and for [8] we locked the affected parameters to zero. As before, we estimated the 3DoF parameters of the camera relative to the ground plane only once for all sequences for our method. The method in [8] requires metrically accurate trajectories. Therefore, we applied the scale factor estimated by our method to the 3D trajectories estimated with ORB-SLAM. The results are similar to the ones from the simulation experiment (Section V-A). Both methods agree on the $x_{y}$ and $x_{\theta}$ parameters, while there are slight discrepancies in the other parameters. Looking at the deviations, both methods are consistent with their estimates. We argue that these differences are due to the time delay estimation in [8]. Overall, our method provides calibration parameters closer to the measured ones.

VI Conclusions

In this paper, we presented a new method for the extrinsic calibration of multiple sensors, suitable for automatic execution on mobile robots. In particular, we first formulated a least-squares problem to estimate the 2D calibration parameters of two coplanar sensors from incremental motions. Next, we relaxed the coplanarity requirement by first estimating the 3DoF transformation relative to the ground plane, and then projecting planar motions into the common plane. Then, we extended the two-sensor case by considering all sensor pairs in a joint least-squares framework. Our formulation allows to accurately estimate the 6DoF calibration parameters of multiple heterogeneous sensors under the assumptions that: i) the ground plane can be observed and ii) accurate per-sensor motions are available. A scale factor is also considered as an estimation parameter during the motion-based calibration and, therefore, our method can also handle monocular cameras. Finally, the proposed approach has been validated with simulated data and assessed in both indoor and outdoor scenarios, supporting its practical application and enhancing the performance of a state-of-the-art motion-based approach. Currently, we are investigating closed form solutions to estimate sensors’ time delays.