{'fraction_non_alphanumeric': 0.02494880546075085, 'fraction_numerical': 0.006262798634812287, 'mean_word_length': 4.614736035259174, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
RGBiD-SLAM for Accurate Real-time Localisation and 3D Mapping

In this paper we present a complete SLAM system for RGB-D cameras, namely RGB-iD SLAM. The presented approach is a dense direct SLAM method with the main characteristic of working with the depth maps in inverse depth parametrisation for the routines of dense alignment or keyframe fusion. The system consists in 2 CPU threads working in parallel, which share the use of the GPU for dense alignment and keyframe fusion routines. The first thread is a front-end operating at frame rate, which processes every incoming frame from the RGB-D sensor to compute the incremental odometry and integrate it in a keyframe which is changed periodically following a covisibility-based strategy. The second thread is a back-end which receives keyframes from the front-end. This thread is in charge of segmenting the keyframes based on their structure, describing them using Bags of Words, trying to find potential loop closures with previous keyframes, and in such case perform pose-graph optimisation for trajectory correction. In addition, our system allows is able to compute the odometry both with unregistered and registered depth maps, allowing to use customised calibrations of the RGB-D sensor. As a consequence in the paper we also propose a detailed calibration pipeline to compute customised calibrations for particular RGB-D cameras. The experiments with our approach in the TUM RGB-D benchmark datasets show results superior in accuracy to the state-of-the-art in many of the sequences. The code has been made available on-line for research purposes https://github.com/dangut/RGBiD-SLAM.


Introduction
In the last years visual SLAM has become one fertile research topic in the fields of computer vision and robotics. Visual SLAM can be addressed either by feature-based or direct feature-less methods. In a previous step feature-based methods extract pixel locations with high gradients and compute the motion by minimisation of the pixelic reprojection error. Feature-less or direct methods skip this step and compute the motion estimate directly from the raw camera input. Regardless which method is used a complete visual SLAM system usually consists of the following modules: • Camera tracking • 3D reconstruction • Place recognition and loop closure when a zone is revisited One of the earliest and most influential works on real time monocular SLAM was presented by Davison [9]. It is a visual feature-based SLAM system, where the camera motion and position of landmarks of the environment are estimated using an Extended Kalman Filter (EKF) from the image displacement of 1 https://github.com/dangut/RGBiD-SLAM the image projections, or features, of such landmarks. Following the success of this work, several contributions introduced new improvements and functionalities like the inverse depth parametrisation of the landmarks [8], a relocalisation module [44], and efficient scheme for data association and outlier rejection [7] or a generalised projection model to support catadioptric sensors [16].
One weakness of this approach is that following a probabilistic filter paradigm, camera tracking and mapping are performed both at frame rate, which limits severally the number of landmarks which could be reconstructed. Klein and Murray [24] addressed this issue with their PTAM (Parallel Tracking and Mapping), where camera tracking and mapping are separated into two different threads. The tracking thread operates at frame rate and is on charge of estimating the camera position given a fixed map; while the mapping thread operates at a lower rate, by building or updating a map by performing bundle adjustment between a set of keyframes selected from the total number of frames. Mur-Artal et al. presented in [28] a visual SLAM approach addressing the weaknesses of PTAM, like handling of landmarks occlusions or mapping scalability, as well as including a novel loop closure and relocalisation method based on Bags of Binary Words using ORB descriptors.
Compared to feature-based methods like the aforementioned direct methods have the advantage of skipping the feature extraction step which performed only an a small fraction of the image pixels and is frequently based on heuristics. However feature-less or direct monocular SLAM methods are more challenging due to the lack of depth measurements from vision sensors and the high computational burden of performing per-pixel operations.
Maybe the earliest and most influencing work in real-time dense monocular method was DTAM (Dense Tracking and Mapping) presented by Newcombe et al. [32]. The authors propose to reconstruct the depth of the scene in keyframes and at the same time estimate the pose of the camera. The problem of estimating a dense depth map from just RGB images is ill-posed and thus authors used an optimisation framework based on variational methods and implemented on GPU. In [13], Engel et al. propose a direct monocular SLAM method which obtains a semi-dense mapping based on variable baseline matching.
Nevertheless, in the last years visual SLAM has been very influenced by two differential events which have rendered direct feature-free approaches more tractable. First the development of a new generation of GPUs and CPUs, which allowed to design algorithms with high parallelisation; and secondly the advent of depth sensors, which are able to provide a dense depth image of the observed environment. The earlier prompted the development of real-time dense visual SLAM methods, while the later enabled SLAM, both feature-based and direct, using RGB and depth channels from a RGB-D camera.
Though RGB-D cameras have the great advantage of providing dense depth information, with respect to RGB cameras they present a more complex calibration. To perform localisation and mapping with RGB, one only needs to calibrate a conventional camera which nowadays is a relatively simple process by using for example Bouguet calibration toolbox [3]. On the other hand, RGB-D sensors require a stereo calibration between the RGB and depth cameras for depth registration. In addition to that, since the depth sensor is really measuring a pixel disparity and then converts it to depth, this relation has to be calibrated also. RGB-D sensors have this functions hard-coded in the hardware, being able depth computation and registration internally. However the use of a fixed factory calibration does not account for imperfections of particular sensors and a customised calibration might be necessary to reach good accuracy in SLAM with RGB-D sensors.
In this paper we present a complete SLAM system for custom calibrated RGB-D sensors. The developed camera tracking module builds upon our dense RGB-D odometry algorithm [17] which presented the novelty of minimising a geometric error parametrised in inverse depth as it fits better the error model of the depth sensor noise. Relying only on the camera tracking for localisation is prone to produce large deformations in the trajectory and the reconstructed map of the environment as the incremental estimate of the odometry drifts. To correct this we incorporate a place recognition and loop closure module which is based in a Bag of Words scheme using ORB de-scriptors [29]. For the 3D reconstruction of the environment we use a keyframe-based approach where spatially close frames are integrated in one single keyframe to reduce the depth sensor noise. To build the 3D model the point cloud of each new keyframes are concatenated to the map point cloud, avoiding the addition of redundant points. Compared to volumetric approaches, 3D mapping on keyframes is more memory efficient as well as allows to better handle the uncertainty in the depth measures during the frame fusion process. In the experiments we show that our method obtains very accurate results in the trajectory estimation compared against a large selection of the most successful state of the art methods in RGB-D mapping.
Additionally, we present a procedure to completely calibrate RGB-D sensor. The main feature of the proposed method compared to most works in the literature is that in order to calibrate the conversion from the measured depth to the real depth we do not rely in any estimate obtained from the depth camera to generate the ground truth depth. The code has been made available in a public repository for research purposes.
The paper is divided as follows. In Section 2 we make a revision of the existing method on RGB-D SLAM and odometry, as well as RGB-D calibration. Section 3 is dedicated to explain some preliminary concepts of the mathematical representation of rigid body motions. In Section 4, we present the models and the pipeline followed for RGB-D camera calibration. Section 5 presents the different image warping approaches that are used in our SLAM approach. Section 6 explains our algorithm for dense frame alignment required to compute the inter-camera motion and Section 7 describes the complete SLAM system with its different modules. In sections 8 and 9 we present the experiments and the final conclusions.

Related Work
Maybe the best example of a feature-based RGB-D SLAM system is the approach proposed by Endres et al. [12]. To estimate the camera motion they first compute an initial seed by RANSAC-based 3D align-ment of sparse features, followed by a refinement step where they only minimise the point cloud alignment error from the ICP algorithm. The method is improved in [11], including an Environment Measurement Model to prune wrong motion estimates which passed undetected in the RANSAC and ICP steps. The system is able to close loops taking a random sample from a set of keyframes and each frame-toframe motion calculation between candidates for loop closure is parallelised, thus the computational cost is close to the frame rate of the camera, between 5 and 15 Hz.
KinectFusion by Newcombe et al. [31] was a ground-breaking contribution in dense RGB-D SLAM. KinectFusion is composed by two different modules, one for camera tracking and one for dense volumetric mapping. For each new frame first the motion is estimated by frame-to-model ICP alignment of depth maps, i.e., current depth map is aligned with the depth map raycasted from a voxelized 3D model. Then, current depth map is integrated in the 3D model using a truncated signed distance function.
The main constraint of KinectFusion is its limitation to small workspaces, which was nevertheless solved in latter works by using a cyclical buffer to shift the volume as the camera explores the environment [1], [41]. However, artefacts are likely to appear in the reconstructed 3D model when an area is revisited. This is solved in recent works by Whelan et al. [42,43] where a loop closure back-end is introduced. This back-end is able to correct potential artefacts in the 3D volume by enforcing the loop constraints not only on the camera trajectory but also on the dense 3D map using a deformation graph.
Dong et al. [10] propose the combination of Kinect-Fusion with a robust odometry algorithm based on keypoint matching and RANSAC alignment between frames. This step provides an initial guess of the sensor pose which is further refined by the ICP step of KinectFusion.
Bylow et al. [5] proposed a method which, as KinectFusion, uses only the depth fusion, but rather than raycasting a depth map from the model for posterior ICP-alignment, the camera is tracked by directly minimising the signed distance function be-tween the current warped depth map and the model surface defined in the voxelised volume. This results in a better accuracy and similar real-time computational performance compared to KinectFusion.
Jaimez and González-Jiménez [21] develop a fast visual odometry method using only depth images. Their method performs frame-to-frame tracking at a camera frame rate of 60Hz and it can be extended to any kind of range sensor.
Following the paradigm of working on 3D models [31], Stuckler and Behnke propose in [36] and [37] converting the RGB and depth images into multiresolution surfel maps by using a voxel octree representation. Each surfel maintains a shape-texture descriptor, which guide data association between surfels in different maps during camera pose estimation. To alleviate the odometry drift they register the current frame with respect to the latest keyframe. A new keyframe is inserted when camera motion w.r.t. last keyframe is large enough. They also propose a loop closure technique where loop closure candidates are randomly sampled from a probability density function which positively weights the selection of spatially closer keyframes.
Kerl et al. [23] propose an RGB-D SLAM system where camera motion is estimated with respect to keyframes, which are switched following an entropybased criteria. They include of a simple but effective loop closure method based on keyframes spatial proximity to further refine the final odometry estimation. The final 3D is obtained by simply joining all the raw depth maps from the keyframes.
Regarding the 3D reconstruction, though RGB-D sensors already provide a dense depth map of the scene, raw depth maps direct from the sensor show a large noise which grows quadratically with the distance of the observed point. For this reason, in order to obtain an accurate dense 3D reconstruction of the environment one needs to fuse multiple consecutive depth maps. As mentioned earlier in this section approaches like KinectFusion, integrate each of the incoming frames in a 3D volume of voxels.
Contrary to volume-based approaches, Meilland and Comport [27], integrate the depth map of each received frame in a set of close keyframes. To track every new frame they use the integrated keyframes to synthesise a reference keyframe from which the new frame is tracked. Though not reporting the use of any loop closure technique their method shows a compelling accuracy in terms of Absolute Trajectory Error, keeping a low trajectory drift.
When comparing both approaches for 3D mapping, volumetric fusion can be viewed as the integration of every incoming frame within a set of concatenated 2D slices, with each slice corresponding to a quantised depth value. This representation in many depth slices can handle occlusions; however it demands large amounts of memory and also, in the process of quantisation and integration of depth measures in the volume, important information about the depth uncertainty is lost.
In keyframe fusion in turn, integration is performed in a single slice represented by the depth map of the frame at which the keyframe was initialised. The main benefits are two. First, a more lightweight and efficient memory representation and secondly, since integration is performed in the frame of the original image space without discretising the depth values, depth uncertainty can be rigorously handled. The main disadvantage of this representation is that, with a single slice, occlusions between elements in the scene cannot be handled. In our keyframe-based approach, since we generate new keyframes periodically this is not a problem. Also it is worth remarking that a keyframe representation comes on handy when performing applying state of the art appearance-based loop closure strategies in visual SLAM, since it provides a direct mapping of the 3D scene to the RGB image of the keyframe. Furthermore, in our keyframe integration scheme, as in our previous work on direct RGB-D odometry, we take advantage of the usage of inverse depth maps instead of depth maps, which allows us to use a constant uncertainty-based tolerance for the integration of a new inverse depth value.
Having a custom calibration for one specific camera is key step to achieve a good precision in SLAM. For this reason another objective of this paper is to obtain a RGB-D SLAM system which can use a custom calibration, as well as providing a calibration pipeline to obtain an accurate calibration. Note that the calibration of an RGB-D sensor is more complex than for RGB ones since it requires first the spatial trans-formation between the RGB and depth cameras, and secondly a function to compute the depth from the raw disparity maps computed by the camera. In the literature, Smisek et al. [35] perform a stereo calibration of RGB and depth cameras, using the infra-red channel of the later. For the correction of a spatial depth distortion they use images of a flat surface, and correct the error wrt. the best fitting plane. Herrera et al. [18], propose a similar approach, but instead they use the depth image for the stereo calibration and perform the distortion correction in the disparity space rather than depth. Di Cicco et al. [6] propose a method for depth distortion correction where some images of a wall are taken. Then a plane is fitted to the pixels in the centre of the images, under the assumption that these pixels suffer less distortion. Finally a pattern of depth multipliers for distortion correction is obtained by regression. Unsupervised approaches for the depth camera intrinsic calibration are proposed by Zhou and Koltun [46] and Teichman et al. [39], where the patterns for depth undistortion are computed taking a ground truth trajectory and map built from a SLAM algorithm. In general, the described approaches tend to use the depth images we want to undistort to generate the training examples to compute the correction parameters for those distorted images. In [45], Zeisl and Pollefeys are aware of this issue and propose an auto-calibration method, where a ground truth of camera poses and sparse map points is generated by a Structure from Motion algorithm using solely RGB images. In this work we follow the indication in [45] of generating the ground truth for depth undistortion from RGB data only. However in our case the calibration process is rather supervised by the user. Though it may be a slower and more complex process, we observe that accuracy is increased, since a supervised method allows first to introduce some measure of the scale such as some distance between points, and secondly to control the scene to be only a textured planar surface so that we can generate a dense depth image ground truth from camera poses computed by photogrammetry and the plane constraint.

Preliminaries: Rigid Body Transformations
This section introduces some concepts of the theory of rigid body transformations which are widely used in this work. For a more complete explanation we refer the reader to [30]. The motion of a rigid body in the 3D space is described by a rotation and a translation with respect to a reference frame. It is usually denoted by a 4 × 4 matrix: which reads as the translation and rotation of reference frame C, which is attached to the body, wrt. another reference frame W .
are respectively the rotation matrix and the translation vector. The rigid body transformations belong then to the product space of R 3 and SO(3) and is denoted as the Special Euclidean Group, The group SO(3) induces an over-parametrisation of the rotations by 3x3 matrices with only 3 DoF, which is a source of problems in optimisation or handling rotation uncertainty. To alleviate this it is usual to switch between a matrix and an axis-angle representation consisting on a 3D vector, θ ∈ R 3 . This representation is also denoted as the Lie Algebra so(3) of SO(3). The exponential map exp : is used to convert a rotation from axis-angle into a matrix representation: where [·] × denotes the skew symmetric 3x3 matrix built from a 3D vector. A rotation matrix is converted to the axis-angle representation through the logarithmic map log : The notion of exponential and logarithmic map can also be defined in the Group SE(3), where an element in its Lie algebra se(3) is given by an screw motion ξ = (v, θ) T , with 6 DoF (3 for rotation θ, and 3 for translation v) : where Q = 1 0 exp(τ θ))dτ and yielding t = Q(θ)v. To avoid the rotation having an effect on the translation, in this paper we redefine the exponential and logarithmic maps to act separately in the rotation and translation, that is: To simplify the notation, from here in advance exponential and logarithmic functions appearing without the group sub-index will refer to the map of the defined above which matches with the function arguments.

Camera model and calibration
A RGB-D sensor is composed by one RGB and one infra-red (IR) camera, as well as a infra-red emitter which projects a dense infra-red 2-dimensional pattern. The projected IR pattern is perceived by the IR camera, and by matching it to a reference pattern obtained at a known depth, a dense disparity map can be computed in a similar way as conventional stereo systems can do sparsely in parts of the image with high intensity gradient. A RGB-D sensor can provide RGB, IR and depth images, where the IR and depth images share the same reference frame and RGB and IR cameras are separated and thus the corresponding images have different coordinate systems. RGB-D cameras include a built-in calibration in their hardware which allows for fast registration of the depth image to the RGB frame. However, this factory calibration does not account for possible imperfections for a particular camera, and thus a custom calibration is necessary to achieve better projection and registration model for a RGB-D sensor.

RGB and IR intrinsics
Both cameras can be modelled by a typical pinhole camera model with the addition of radial and tangential distortion correction [4]. A scene point X = (X, Y, Z) T is first projected on a normalised plane obtaining the undistorted coordinates mu = (X/Z, Y /Z, 1) T . Then the distortion correction is applied: The vector k = (k1, k2, k3, k4, k5) T encapsulates the distortion parameters. This step is followed by a linear projection from the normalised plane in the final image plane: where K is the pinhole model calibration matrix including the focal length f = (fx, fy) T and the principal point c = (cx, cy, 1) T . All the previous steps can be encapsulated in a projectivity function p = π(X).
Image points can be also lifted to the 3D world with the inverse model: Note that in this last equation, for X arbitrarily Z = 1 since depth cannot be recovered from one single projection. It is necessary to know the depth Z(p) in a pixel to scale X and obtain a full 3D point.
The depth camera provides a depth image which is obtained from the disparity map using a generic conversion function. Thus, additional calibration parameters are needed to improve the depth map accuracy for a specific camera. We apply the following model to the inverted raw depth map Wm = 1 Zm measured by the camera : where β0 and β1 are the calibration parameters, and p0 = (4, 4, 0) T is a pixel shift caused by the correlation window used internally by the depth sensor to compute the disparity [25]. The purpose of this function is to correct the the disparity-to-depth conversion computed internally by the sensor using factory constants. However, the depth sensor introduces also a spatially dependent error. To correct this distortion we apply an undistortion function as follows: where D 0 (p) and D 1 (p) are spatially varying offset and multiplying factor respectively. Each of these terms is calibrated with 8 parameters which are the coefficients of a polynomial in the 2D image coordinates: Di(p) = qi,0 + qi,1r 2 + qi,2r 4 + qi,3r 6 + qi,4mx + qi,5my+ where m = K −1 p and r = m 2 x + m 2 y .

RGB and depth cameras extrinsics
The relative pose between the color and the depth camera frames is denoted by the spatial transforma- It is required to register the inverse depth image in the color camera frame.

Calibration pipeline
The calibration of the RGB-D camera is performed in two steps. The first step consists in a conventional calibration of the IR and RGB cameras intrinsics and stereo transformation using the Bouguet calibration Toolbox [3]. With this procedure we obtain the RGB and depth intrinsics f {C,D} , c {C,D} , k {C,D} and the spatial transformation T D C between the RGB and depth camera frames.
The second step in the calibration is the computation of the depth intrinsics for the correction of the raw inverse depth maps. For this step, we need M synchronised depth and RGB shots taken at different distances from a wall filling the whole image with a checker-board pattern and some posters attached to it. The real inverse depth map is computed by photogrammetry from the RGB images. After corner extraction, image matching and applying bundle adjustment techniques we obtain the pose T W C i of each of the M RGB frames and a 3D sparse point cloud. Since the dimensions of the checker-board pattern are known, the estimated variables are in the true scale. The plane Π = (nW , d) which best fits the point cloud is also computed.
With the poses T W C i , and the stereo camera intrinsics from the previous step all the pixels in the inverse depth image can be lifted and represented in the global reference frame: And by applying the restriction of all the image points being contained in a plane n T we can compute a ground truth for the inverse depth images: Then, the intrinsic depth correction parameters are computed in two steps. In the first step we compute only the parameters of the linear relation: In the second step after applying (12) we compute the coefficients of the distortion polynomials: The reason of optimising in two separate steps is obtaining a calibration where we can choose whether use or not use the spatial distortion correction.

Warping functions
A key element in direct methods for visual slam is the image warping functions. These are need for distortion correction or to obtain representations of the scene in some image as viewed from another different pose. Depending on what we need to achieve we use different warping approaches.

Inverse warping
Given a source image IA and a destination image IB, inverse warping consists into mapping every pixel p ∈ Ω in the destination image to a position in the source image, and then obtain the corresponding value by bilinear interpolation at the source: where the fw is the warping function. We use inverse warping for undistortion of intensity and inverse depth maps, where fw(p) = Kh(K −1 p).

Forward registration
Forward registration is required when we need to represent an inverse depth map in another reference frame, like for example the one of the RGB camera. Generically, given the spatial transformation T A B and assuming that images have been undistorted, a pixel p in a depth map WA at source frame A is mapped to a pixel pB in a destination frame B through the following equation: whereR B A naive approach would be to compute pB, round to the nearest pixel location and then update the inverse depth value at this location. Unfortunately this approach is liable to produce holes in the new registered depth map. Instead we use a solution based on the relief texture mapping algorithm [33], performing the registration in two steps.
In the first step the depth image is registered to an intermediate frame B t applying only the shift due to translation: From this equation we obtain: to prevent the arising of holes we map the whole pixel surface to the frame B t . Adding δp = (δpx, δpy, 0) to p in (21) and assuming that WA(p) = WA(p + δp) ∀ ||δp|| ∞ ≤ 0.5, we get: Taking the four pixel corners at frame B at δpx ± 0.5 and δpy ± 0.5 and rounding the corresponding pB t + δpB t to the nearest integer, we obtain a square containing all the pixels where wB t should be mapped at frame B t . For every pixel ps in this square WB t (ps) is updated if wB t is the greatest or the first inverse depth value computed in that pixel.
In the second step, the remaining rotation is applied to the warped frame. Noting that theR B A represents an homography, and thus a linear mapping, inverse warping can be applied taking pixels in the destination frame B:

Inverse geometric warping
If inverse depth and intensity maps are available at both frames A and B, inverse warping can be performed by using (20) as the warping function. Note that since we are performing inverse warping, the roles of the frames A and B in this equation is interchanged (A is the destination, and B is the source). This is computationally faster than forward warping and is useful, for example, to fuse inverse depth maps taken at different poses. The drawback is that occlusions cannot be handled and thus the motion between frames must be small. Given that the pixel position pB computed by (20), the images at frame B * warped towards frame A are computed by the following relations:

GPU implementation
All the warping functions have been implemented for CUDA capable NVIDIA GPUs. In the inverse warping functions, the mapping of source images to the GPU texture memory allows for efficient interpolation via texture fetching. In the case of the translation-only warping step in the forward registration procedure, a pixel in the destination images might be accessed at the same time by two or more GPU threads for a writing operation. In order to avoid race conditions and serialise the access for writing we use atomic operations.

Dense motion estimation
Dense frame alignment consists into estimating the camera motionT A A ) between two frames A and B by pixel-wise minimisation of the photometric and inverse depth error between both frames.

Photometric and geometric constraints
Let us denote two camera frames as A and B, at instants t and t + ∆t respectively. Given the intensity images IA and IB, and inverse depth maps WA and WB defined over the image domain Ω ⊂ P 2 , for an image point p = (u v 1) T ∈ Ω in frame A, the following constraints hold: where ∆Xp and ∆p are respectively the scene and the optical flow of one point from frame A to B, and e T z = (0 0 1). The constraint in intensity assumes constant illumination of one scene point. The second constraint models the change of one point's inverse depth due to the relative motion along the optical axis of the camera.
It can be verified that the scene and the optical flow satisfy the following equation: With this relation and computing the first order Taylor expansion with respect to the scene flow we get the linearised constraints: Since we are estimating the motion between frames we are working under the assumption of a rigid scene. Assuming a small motion described by the rotation and translation pair ( where [·] × denotes the antisymmetric matrix from a vector and

Iterative optimisation
The inter-frame motion is computed iteratively in a coarse-to-fine manner using image pyramids, performing a given number of iterations at each pyramid level. At the start and every time we step down in the image pyramid, we downsample {IA, WA} and compute the gradients {∇IA, ∇WA} at current pyramid level. The initial estimate of the rigid motion is set to the identity,T A B (0) = I, unless an initial guess is provided.
At each iteration γ, first intensity and depth maps in frame B, {IB, WB}, are warped towards frame A taking the up-to-date motion estimateT A B (γ) . By using the reverse warping approach described in Section 5.3 we obtain the warped images {I (γ) These images are downsampled to the pyramid level at which current iteration is taking place. The formulae for the residuals at current iteration are given by: The next step is the computation of the update of the motion estimate between frames. This is achieved by optimisation with Iteratively Reweighted Least Squares algorithm (IRLS) [19], which results in a linear least-squares problem to be solved at each iteration: The weighting function ω(x) depends on the assumed probability distribution of the residuals. µ {I,W} and σ {I,W} are the location and scaling parameters which capture the bias and the uncertainty in intensity and inverse depth residuals, and allow for normalisation of residuals in different magnitudes. The details on the computation of the weighting function and the location and scaling parameters are presented in the next subsection 6.3. To account for noisier depth measurements from surfaces more oblique to the camera, λn(p) = n T p ||p|| is used to down-weight the inverse depth residuals, where the normal n T is computed as: where α is the normalising constant. The minimisation of (36) leads us to solve the following linear system: and after computing ξ A B (γ) we update the interframe motion estimate for the next iteration: The covariance of the motion estimate is computed as the inverse of the Hessian. The Hessian depends on the intensity and depth maps gradients computed in the frame A. Due to sensor noise it is possible that some pixels provide misleading information making the Hessian larger, which would yield a wrongly optimistic covariance matrix. In order to prevent this, once we have computed the final iteration for the motion estimate, we estimate the Hessian by performing an extra iteration where the motion estimate is not updated and where {IA, WA} are applied a bilateral filter prior to the computation of their gradients:

Error distribution and robust optimisation
The underlying assumption behind the IRLS method is that the residuals we want to minimise are modelled by a probability density function (pdf). This function is related to the weights ω(x) we apply in the cost function. The purpose of IRLS is that, by selecting a pdf with heavier tails than the typically assumed Gaussian distribution, gross errors can be accommodated and and thus we gain robustness to outliers. In this work we use the student t-distribution whose pdf is given by: with x i = ri−µ σ , and where µ, σ are the location, and the scale of the pdf, while ν is a characteristic parameter of the Student's t-distribution and denotes the degrees of freedom. In robust estimation the parameter ν is associated to the achieved robustness. Taking the logarithm of (41) and eliminating constants we obtain the log-likelihood function: The location and scaling parameters of the residuals are obtained by maximisation of their joint loglikelihood function. That is: Taking the derivatives with respect to µ and σ we obtain respectively the following equations: where ω(x, ν) = 1 x ∂ρ(x, ν) ∂x . This system of equations can be solved iteratively using the previous estimates of µ and σ to compute ω(xi, ν) before each iteration.
Note that up to this point we have ignored the estimation of the parameter ν. Typically this parameter is set to ν = 5 which provides enough robustness and also behaves well in case the residuals have a Gaussian distribution. However, ν can be estimated as well, which allows to adapt the required robustness to the distribution of current residuals. By maximising (42) with respect to ν [26] we obtain the following non-linear equation: where φ(y) denotes the digamma function. Note that ideally robustness parameters should be estimated at the same time as location and scaling. However, unlike location and scaling equations, fixing the weights in (46) with the previous estimate does not produce a linear system, requiring from some sub-iterations to compute ν within current iteration. Moreover fixing ν in the weights to its previous value produces a slow convergence. However letting ν to be an unknown in the weights implies to perform a sub-loop of GPU reduction operations within the loop in which the distribution parameters are estimated. Therefore, to keep computational cost low we estimate first the final location and scaling values for intensity and depth residuals fixing ν = 5, and then we solve (46) with the computed values using the bisection method assuming ν ∈ [2, 10] with a lower ν implying more robustness to outliers. This is performed both for the intensity and inverse depth residuals yielding νI and νW , In dense RGB-D tracking the main source of outliers are the high inverse depth residuals caused by occlusions which are due to dynamic objects or areas which simply become revealed or hidden by the camera motion, i.e., outliers are caused by large inverse depth residuals rather than intensity ones. For this reason and because too much robustness can lead to a bad or slow convergence to the solution, we select the highest of the computed νW and νI for the intensity residuals. νI = max(νI, νW ). (47) A naive computation of the location, scaling and robustness parameters would use all the residuals. Given that every residual in intensity and inverse depth corresponds to a single pixel, the number of samples would be extremely large, more than 300000 samples at the highest resolution level, which would involve a high computational cost. To reduce this burden we compute the scaling parameters in a sample with a maximum size of N = 19200 points obtained by systematic selection. This sample size guarantees first an integer stride for pixel sampling at all the considered resolutions; and secondly, after leading an statistical analysis, it also guarantees a relative precision of 0.03 with a confidence level of 99.7%.

Dense frame covisibility ratio
Computing the covisibility score between two given frames is a key element to prune redundant frames, selection of reference frames or keyframes, or computing visibility graphs. In this work, taking advantage of the availability of dense depth maps, we use a covisibility computed from all the pixels in the image. Given two frames A and B and the camera motion estimate between them ( R B A , t A B ) ∈ SE(3), we transfer pixels from A to B using (20) and compute the following pixel subsets: SH is the subset of pixels in frame A with valid inverse depth values (no-hole pixels) and SV the set of visible pixels when transferred to frame B. A pixel is tagged visible if it is within the image bounds, and it is not occluded in frame B. A pixel is considered occluded if its mapped inverse depth value is above an uncertainty allowed by the sensor noise model. Then the covisibility ratio is computed: This procedure is repeated switching the role of A and B, and then we select the minimum ratio as the final covisibility ratio.

RGB-ID SLAM system
Our RGB-ID SLAM system is implemented in two CPU threads running concurrently. One thread executes a front-end for every incoming frame from the RGB-D sensor, performing camera tracking and fusion of inverse depth measurements in a single keyframe. For both of these tasks the CPU thread calls functions which execute in GPU where pixelwise operations can be easily parallelised. The second thread executes a back-end in charge of managing the keyframes passed by the front-end through a buffer. First it performs a segmentation of the keyframe in juts using the map of normals. Secondly for every keyframe it obtains different BoW histograms and stores the keyframe in a database. Finally it attempts to close loops comparing the last processed keyframe against the keyframes in the database. All the tasks of the back-end thread are performed in CPU except for a dense alignment in GPU between keyframes when a successful loop is detected. In order to allow for simultaneous access to GPU by both threads without blocking, each thread calls its GPU functions on its own CUDA stream which is executed asynchronously.

Camera tracking
For every upcoming frame {I k , W k } we have to compute the rigid body transformationT rf k which best aligns it with a given reference frame {I rf , W rf }. This motion is computed by the dense RGB-D alignment method described in Section 6. We provide an initial guess of the motion estimate using a constant velocity model, i.e.,T rf Once we compute the motion estimate between the reference and current frameT rf k and its covariance Σ rf →k rf , we compute the sequential odometry constraint as it will be required in the latter pose-graph optimisation step.
where the computation of the Jacobians of the transformations is detailed in the Appendix A.

Keyframe fusion
After the visual odometry, the inverse depth map of current frame is fused in the last selected keyframe. The criteria for keyframe selection is the same covisibility criteria as used for the reference frames taken for odometry. We must note that generally keyframes do not have to coincide with the odometry reference frames, since we allow for using a different covisibility threshold. Normally it will be set to a lower value for keyframe switching, i.e., we will initialise new keyframes at a lower rate than reference frames, using a visibility threshold of 0.7. Once a keyframe is initialised, the inverse depth maps of the incoming frames are fused with the keyframe's. We not only fuse the incoming frames but also the frames recorded prior to the keyframe selection in order to integrate inverse depth measurements in parts of the image no longer observed after keyframe initialisation. To do so, frames are stored temporally in a buffer and erased once integrated in the later keyframe. Every time we integrate an incoming frame, we integrate also the frame closest in time of the ones remaining in the buffer. Fusion in keyframes is performed in GPU using the reverse warping approach detailed in section 5.3, where the frame to be integrated is warped towards the keyframe. Note that since we are using a reverse warping scheme parts of the scene with missing depth measurements in the frame from which the keyframe is initialised cannot be reconstructed from subsequent frames. We find this a minor disadvantage given first, that reverse warping in CUDA is computationally cheap and easy to implement and secondly, that the generation of many keyframes generally guarantees that zones missing in one keyframe can be observed in another one.
Once the frame has been warped we have to integrate every new inverse depth value for every pixel in the keyframe. To do this we first verify that the pixels on both the keyframe and the warped frame correspond to the same scene point by performing the same covisibility check as in Section 6.4. If a pixel passes this test we update its inverse depth value W kf (p) as well as the corresponding weight C kf (p) in the keyframe as follows: A p 2 stems from the propagation of the inverse depth variance through (20). When a new keyframe is initialised the values of the weight map C kf is reset to 1 and the last keyframe is stored in a buffer, awaiting to be processed by the back-end in the second CPU thread.

Keypoint extraction and BoW histograms
Loop detection primarily relies on the DBoW2 place recognition system based in Bags of Binary Words proposed by Gálvez-López and Tardós [15]. Mur-Artal et al. [29] used this approach with ORB descriptors obtaining a quite reliable and fast loop closing system. As the authors suggest, we compute a pyramid of 8 levels and a relative scale of 1.2 between levels from the intensity image, and then we extract FAST keypoints at each pyramid level. In order to distribute the keypoints uniformly over the image we divide the image in cells and set the firing threshold for the FAST extractor adaptively to extract the desired number of points at each cell. Once the keypoints are extracted, they are lifted to 3D points using the inverse depth map of the keyframe and ORB descriptors are computed in the intensity image. For the classification of the ORB descriptors in a BoW we use the same dictionary as in [29]. For each keyframe a BoW histogram is computed with the frequency of each word in the given image. Finally the processed keyframe is stored in a keyframe database from which keyframes are queried in the loop closure process.

Superjuts and entropy of normals
For every keyframe which is passed to the back-end we generate a map of the degree of disarrangement of the elements in the scene,based on the entropy of the normals. For a given keyframe, we generate a normal map from its inverse depth map. Then we perform an object-like segmentation based on the algorithm for RGB images developed by Felzenszwalb and Huttenlocher [14], and which is similar to the methods presented in [40] and [22]. The main difference is that we perform the segmentation directly on the keyframe pixels rather than on a lifted point cloud.
A graph G = (V, E) is constructed, where each vertex in V represents a 3D point X i lifted from a pixel location and the edges E represent the neighbouring relations between points. Edges are inserted between points with adjacent pixel locations on the image. Having the normals at every point, a weight w ij for the edge joining vertices i and j is computed: where the squared weight is applied to convex edges, reflecting the fact that convex regions usually contain points belonging to the same object and concave regions are likely to arise in frontiers between objects. After computing the weights, the segmentation algorithm is run and essentially groups points sharing edges with low weights in the same superjut. A parameter k must be tuned such that the higher k, the larger segments will be obtained. We set this value always to 0.6. After the segmentation we compute a histogram of normals F (n). Each of the M bins in the histogram corresponds to a Voronoi cell associated to one point out of M points uniformly distributed over the sphere. Since there is no analytical solution to the problem of evenly distributing M points on the sphere for any M , we use a simple approximate solution based on the golden section spiral [2,34], which results in a discretisation with bins covering areas of similar size. Finally given the histogram of normals of a jut, we can compute its entropy: (56)

Loop closure
After extracting keypoints and computing the BoW histogram the keyframe database is searched for potential keyframe candidates for loop closures. Given two keyframes, DBoW2 returns a similarity score based on the distance between their BoW histograms, which we normalise with the score between the currently querying keyframe and the previous keyframe. If this score is above a threshold for any of the histogram comparisons a potential loop candidate is stored. In order to ensure that potential loop closures are obtained between temporally distant keyframes we establish a minimum keyframe separation. For every possible keyframe candidate we match its 3D points with the 3D points of current keyframe using their ORB descriptors. Then we perform a geometric validation of the detected loop by computing a rigid motion estimate applying the method described in [20] for alignment of 3D point clouds in a 3-point RANSAC scheme. If more than 10 points agree with the computed motion between keyframes, and the convex hull spanned by these points takes up more than 5% of the image, the loop is tagged as correct. In order to compute the loop constraint, containing both a motion estimate and its covariance, we align both keyframes using the approach described Sec. 6 where we pass the RANSAC rigid motion estimate as initial guess.

Pose-graph optimisation
When a new loop is detected it is added to a cluster containing loop constraints which are tem-porally close, allowing a maximum separation of 10 keyframes between any two of the keyframes in the cluster. Pose-graph optimisation is applied when one of the clusters is no longer accepting new loop constraints, updating the graph with the new constraints in the cluster, as well as the poses and odometry constraints computed by the frontend, since the last graph update. Then, given a graph with a set of odometry constraints, S = { (1, 2), ..., (N − 1, N )}, and a set of loop constraints, R = {(i1, j1), ..., (iL, jL)}, we want to find the tra- which minimises the following cost function: where Ω ij is the information matrix of the constraint and the error term eij(x) represents the residual between a constraintT i j and the relative motion between the trajectory poses T W i and T W j it is connecting, expressed with a minimal parametrisation: From this error term, the Jacobians and information matrix required for each constraints can be computed as detailed in Appendix B.
To speed-up the process we perform the optimisation in a multi-layered scheme. In the first layer, odometry constraints between consecutive frames, are substituted by constraints between consecutive keyframes, which are also computed by the tracking front-end, and thus we only optimise the poses of the keyframes. In the second level the keyframe poses are fixed, and the poses of the rest of the frames are optimised by enforcing the odometry constraints.

Experiments
We first compare our method quantitatively with the state-of-the-art in terms of trajectory estimation. Then we perform a qualitative evaluation of the 3D reconstructions both in the Freiburg and our own RGB-D sequences, and finally we measure the computational performance, both of the front-end, in charge of camera tracking, and the back-end, performing keyframe processing and loop closing.

Trajectory Estimation
We first compare our method to state-of-the-art visual odometry and SLAM approaches with RGB-D systems. For the comparison we use the TUM benchmarking dataset [38]. The evaluation has been carried out taking the error metric of the Absolute Trajectory Error (ATE) in meters (see Table 1). However it is worth noticing that there exists a large improvement in the f r2/desk dataset. Apparently the reason is a correction we applied to the depth maps of the Freiburg 2 sequences by a constant scaling factor. This correction was in principle reported as done by the authors of the dataset. However it seems from the observation made by Mur-Artal et al. in [29] of a scale bias in the RGB-D SLAM system of Endres et al. [11], that it has not been applied in the Freiburg 2 datasets available for download.
The table shows that our method outperforms most of the odometry and SLAM methods in the literature and it is close to the ElasticFusion of Whelan et al. [43], which shows the best performance in the state-of-the-art in the datasets f r3/of f ice and f r3/nst. Among the methods performing only odometry, ElasticFusion in tracking-only mode outperforms them and even some of the methods performing loop closure. However in other datasets our method shows overall the best performance both in trackingonly mode and with loop closures. It is worth noticing that challenging datasets like f r1/desk2 and f r1/room, showing a fast camera motion, our method provides a great precision.

3D reconstruction
In this section we perform a qualitative evaluation of the 3D models reconstructed of our SLAM system. To build the 3D model of the scene we concatenate the point clouds of different keyframes in real-time. To avoid unnecessary redundant points we only concatenate partial point clouds which correspond to novel parts of the scene between consecutive keyframes. Redundant points due too loop closure are eliminated at the end of the execution of our method by performing a voxel grid filtering with a voxel size of 1 cm. In Fig. 2 we show the 3D models of the evaluated TUM datasets, showing a great detail of the reconstructions.

Computational performance
The computational performance of our system has been evaluated in a laptop PC with an Intel Core i7-4710HQ CPU at 2.50GHz, 16GB of RAM and a nVidia GeForce GTX 850M GPU with 4GB of memory. As shown in Table 2, the execution time of the front-end varies depending on whether the system operates in tracking only-mode, i.e., with the back-end disabled or with its complete functionalities. This is due to the division of the GPU resources upon successful loop closure between the visual odometry and integration modules from the front-end, and the keyframe alignment module between loop keyframes from the back-end.
In Fig. 3 we show graphically the computational cost of the back-end pipeline, which has been dissected in the different processes involved. As we can observe there is a nearly constant cost for keyframe segmentation and descriptor extraction, while the cost for loop closure detection and pose-graph optimisation shows a great variability. The reason is that most of the cost of loop closure detection is governed by the step of geometric verification and loop constraint computation, which is run only on a reduced list of loop candidates delivered by the efficient appearance-based BoW algorithm for loop detection. Red dashed line represents the time it takes the front-end to deliver a new keyframe to the back-end, assuming that it is able to process al the frames in the sequence at the frame rate of 30 Hz. We can observe that in datasets where the camera is moved slowly in a loop around the same scene, the back-end processing time is widely below the average keyframe rate. However in datasets like f r1/desk2 or f r1/room, where the camera moves quickly and/or the area being mapped changes a lot, the back-end Table 1: Absolute trajectory error (RMSE, median and max) in meters of our method with and without loop closure and comparison with state-of-the-art approaches. For a given error measure in a dataset, we show in bold the best approach in the state of the art and ours if it is the absolute best.  [5] 0.035 [43] 0.022 struggles to operate at keyframe rate. These observations are embodied numerically in Table 3.

Sensor Calibration and Depth Correction
In this section we evaluate the calibration model and the procedure proposed in Section 4 to compute a custom calibration for a specific sensor. The calibration is done in two steps. In the first step we take pairs of RGB and IR images compute the RGB and IR intrinsics as well as the rigid transformation between the RGB and depth reference frames, using the Bouguet Toolbox [3]. In the second step we calibrate the disparity-to-depth conversion of the depth sensor.
To do this, we take 18 shots of synchronised RGB and depth images with the camera pointing a textured wall. From the RGB images we have generated a ground truth, by manually selecting and matching some dozens of points across images, and then performing a Bundle Adjustment with the Photomodeler software to obtain the camera poses as well as a 3D point cloud from which the best fitting plane is computed. Then, using (16) with the calibration parameters computed in the first step, we generate a dense inverse depth map which is used as ground truth in the optimisation function used for calibration of the disparity-to-depth conversion. To test the accuracy of our method we have selected the odd images for calibration and the even ones for validation. The comparison of our calibration, with and without   In the datasets where the average acquisition time is not shown, it means that it is out of the plot bounds (above 600 ms). This usually occurs in datasets when the camera executes a loop around some scene, which normally involves that keyframes are switched with low frequency.
depth spatial distortion correction, and the factory calibration is shown in Fig. 4. We can observe in the figure, that there exist a bias in the raw depth map which is completely corrected by calibrating only the linear model. However, after this adjustment it still remains a zero-mean error with a spatial pattern, which is corrected by the calibrated model of the depth spatial distortion. Note that with a custom calibration the depth error is reduced in both the calibration and the validation images, which indicates that the disparity-to-depth conversion is well modelled by our calibration model. To evaluate the performance of our RGB-iD SLAM system with a custom calibration we have run our method on one sequence taken with synchronised RGB and depth streams but with unregistered depth maps. The obtained reconstruction is shown in Fig.  5

Conclusions
In this paper we have presented a complete RGB-ID SLAM system which combines an accurate dense and direct RGB-D odometry system and state-of-the-art loop closing procedure to achieve a great accuracy in the estimation of the trajectory. The algorithm also performs the integration of the tracked frames into keyframes with smooth inverse depth maps. The fu-   sion of depth measurements in keyframes also allows to achieve a greater accuracy in the computation of relative camera transformation at loop closure. The concatenation of the novel parts between consecutive keyframes allow to obtain precise 3D models of the environment. In the comparison with some of the best state of-the-art methods for RGB-D mapping our method shows to have results with greater or similar accuracy in most of the tested datasets. Our method is able to work in real time with a mid-range laptop GPU. Our system is also to work with custom calibrated RGB-D sensors without need of hardware registering the depth image to the RGB frame. With this in mind we have also proposed a pipeline for accurate calibration of RGB-D cameras. The code for the RGBiD-SLAM system has been made available publicly 2 .

A Covariance propagation in composition of spatial transformation
Let us define the following composition operation of spatial transformations T respectively expressed in the reference W , and we want to compute the covariance Σ AB A . To compute the Jacobians we must apply differential perturbations in the side in which the covariance is referenced. In this case all of them are left perturbations. Taking separately the translation and rotation we have: Applying that exp(δθ) ≈ I + [δθ]× and [a]×b = −[b]×a in (59); and R exp(θ)R T = exp(Rθ) and log(exp(δθ1) exp(δθ2)) ≈ δθ1 + δθ2 in (60) we obtain the following expression for the covariance propagation: B Pose-graph Jacobians and Information matrix Frequently in optimisation problems, the measurement contribution to the error term is expressed through an addition or subtraction, and thus the information matrix of the constraint is trivially the inverse of the measurement's covariance. However in a pose-graph problem the error is expressed as: and thus we need to compute the derivative of eij wrt.T i j . Assuming that the covariance ofT i j , Σ ij i , is computed on the left side, we will compute the derivative applying a perturbation at the left ofT i j : The derivative of the error wrt. the relative pose measurement remains: where ∆Qij = 1 0 exp(τ log(∆Rij))dτ . And hence we obtain a expression for information update: To obtain the Jacobians wrt. to the state of the optimised nodes we proceed similarly. The incremental updates are to be applied in the form of transformation matrices by post-multiplication. Thus in this case the perturbations to compute the Jacobians are applied on the right side of the absolute camera poses in (64). Thus we obtain: