{'fraction_non_alphanumeric': 0.02848799943024001, 'fraction_numerical': 0.009044939819101204, 'mean_word_length': 4.787223080886141, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Elasticity Meets Continuous-Time: Map-Centric Dense 3D LiDAR SLAM

Map-centric SLAM utilizes elasticity as a means of loop closure. This approach reduces the cost of loop closure while still provides large-scale fusion-based dense maps, when compared to the trajectory-centric SLAM approaches. In this paper, we present a novel framework for 3D LiDAR-based map-centric SLAM. Having the advantages of a map-centric approach, our method exhibits new features to overcome the shortcomings of existing systems, associated with multi-modal sensor fusion and LiDAR motion distortion. This is accomplished through the use of a local Continuous-Time (CT) trajectory representation. Also, our surface resolution preservative matching algorithm and Wishart-based surfel fusion model enables non-redundant yet dense mapping. Furthermore, we present a robust metric loop closure model to make the approach stable regardless of where the loop closure occurs. Finally, we demonstrate our approach through both simulation and real data experiments using multiple sensor payload configurations and environments to illustrate its utility and robustness.


I. INTRODUCTION
M AP-centric SLAM solutions [1], [2], which have demonstrated their accuracy and effectiveness by fusion-based mapping and deformation-based loop closure, provide an alternative solution to the dominant trajectorycentric LiDAR SLAM [3]- [5]. Due to the benefit that its complexity does not increase over time, the high quality reconstruction can be easily achieved by having redundant scans and fusing them. Despite their success, previous map-centric approaches have a particular set of drawbacks [1], [2], [6]- [8]. Firstly, the original map-centric work is dedicated to a pinhole camera model such as RGB-D [1] or multi-beam LiDARs [6]. Secondly, asynchronous sensor fusion and handling high motion distortion of LiDAR [9] are not well studied within the map-centric framework. This is far from ideal as LiDARbased systems often need to be tightly fused with multi-modal sensors to handle the motion distortion effects. Lastly, the map-centric approach is extremely susceptible to loop closure failure. The sensor data history cannot be stored due to the extensive input data size which implies that, once incorrectly fused, the map is not recoverable.
In this paper, we seek an algorithm that exhibits the desirable characteristics of the map-centric approach, but is able to handle these problematic situations described above. We propose to do this by introducing the Continuous-Time framework into the map-centric approach. The detailed contributions of this paper are as follows: • We introduce a new constraint for local Continuous-Time trajectory optimization to operate like a map-centric approach. • We detail the realization of a new Continuous-Time trajectory optimization approach which is dedicated to the map-centric approach, improving local trajectory estimation accuracy. • We derive a novel probabilistic surfel representation using the Wishart model. To the best of our knowledge, this is the first approach to introduce the Wishart model for surfel fusion. • We introduce a surface resolution preservative surfel matching method for non-pinhole type sensor data. • We implement a LiDAR only sequential metric localization to reduce failure on loop closure. • We demonstrate the solution on various environments and sensor configurations using simulation and real datasets.
Our previous work in [2], [10], [11] presented an approach for LiDAR mapping using probabilistic surfel fusion, a method for map deformation, and a method for robust Visual-LiDAR metric localization respectively. These approaches deal with the issues separately and independently. In this work, we reformulate the map-centric 3D dense mapping problem with our new trajectory optimization method and intensively revise the dense surfel fusion model. Furthermore, in the metric localization stage we weaken the dependency of the vision to reduce the uncertainty from unobservable parameters related to the vision while maintaining the accuracy. Finally, we evaluate, analyze, and demonstrate this new approach in a comprehensive manner across various sensing platforms and across different environments (Figure 1) using simulation and real datasets to reveal its utility and robustness. arXiv:2008.02274v1 [cs.RO] 5 Aug 2020

A. Map-centric Approach
ElasticFusion [1] proposed a map-centric approach for RGB-D cameras that removes pose graph optimization yet performs globally consistent mapping by giving elastic property to the map and directly deforming the map. The concept of map-centric SLAM eliminated the need of a pose graph for globally consistent mapping and converted the time dependency of the global optimization to a space-dependent problem. Also, by confining the tracking and fusion within recent map elements, they drastically reduced the processing time per input frame.
Despite these improvements in ElasticFusion method, some features of their approach are limited to RGB-D sensors and are not applicable to a LiDAR sensor model [10]. Firstly, the original frame-to-frame motion model [1] is not an effective framework for handling asynchronous estimation and severe motion distortion of LiDAR. Secondly, the map fusion method which is based on a pin-hole camera model can be applied only to the LiDAR configuration with low vertical field of view [6]. For example, the pin-hole camera model cannot be applied for the spinning LiDAR system [3], [12] where the scan range often covers 360 degree both vertically and horizontally. Thirdly, conventional LiDAR loop closure models [13] are not suitable for the map-centric approach. The second and third points will be elaborated in further detail in the following sections.
In a broad sense, sub-map dividing and realignment approaches [8], [14] are also similar to the map-centric method due to its property that input frames are fused into the local submaps. While the approaches in this category successfully reduced global optimization cost, the map fusion policy between graph nodes cannot be defined which brings the discontinuity issue in the map maintenance.
Recently, Behley et al. [6] presented an approach analogous to the concept of the map-centric where the new estimations are fused into the global map. Although, the density of their map representation is dedicated for localization rather than dense mapping, their result present an important direction toward the fusible LiDAR map representation along with our previous research on LiDAR surfel fusion [10].

1) Map Representation and Fusion Feasibility:
The fusion between the global map and input frame is an essential step for high-quality map reconstruction and frame-to-model localization. A number of techniques exist for map representation and fusion in SLAM as shown in Figure 2 (a).
Most common dense map representation is using raw point cloud [?]. While it is most simple and efficient way to visualize  the reconstruction result, the point representation is not proper for map element wise fusion and also requires relatively far more number of map elements to make a dense visualization. The approaches that digitize the space [15]- [18] successfully modeled probabilistic fusion of space cells for the free and occupied space categorization by ray tracing. However, the space digitization create non-smooth reconstruction which limits its applications to navigation and obstacle avoidance where map resolution is less important [19]. A mesh map representation provides a smooth object representation but the map update rule could be inefficient as the connections among the vertices are needed to be properly updated [20].
In [21], [22], the authors presented a multi-resolution surfel map representation and fusion for a spinning 2D laser scanner. Although this representation is highly suitable for fast and robust map registration, large size surfels and its z direction volume does not make it ideal for a dense map representation.
On the other hand, a dense surfel map representation [1], [6], [23]- [26], originally designed to render 3D point clouds without a complicated mesh extraction step [27], is ideal for fusion, dense object representation, and rendering as shown in Figure 2 (b). Thus, a surfel representation has been popular in RGB-D reconstruction problems where a rendered synthetic 2D image from the densely reconstructed object is immediately used for other purposes such as localization [28] or dynamic object handling [24]. Despite its advantage, its introduction to LiDAR-based system has not been actively studied in the LiDAR SLAM community due to the difficulties in the data association which will be discussed in the following section.
2) Map Element Matching: Data association is one of the most important components in map fusion. The simplest way is to find a surfel with the closest distance either in Euclidean space [29] or Mahalanobis distance [30], [31]. However, generally it is challenging to control the map surface resolution without discretizing the environment [16].
The projective data association [24] provides an efficient map element matching but it is restricted to a projective sensor model. Behley et al. [6] applied projective data association and frame-to-model registration for the multi-beam LiDAR sensor. Their approach demonstrated that the projective data association can be successfully adapted to multi-beam LiDAR.
While approaches in this category are computationally more efficient, the frame-to-model registration method is not suitable for handling severe motion distortion and asynchronous multi-modal sensor payload (Visual-LiDAR-Inertial). Furthermore, the projective data association is limited to the LiDARs with the projective view model whereas the modern LiDAR systems [32] often cover 360 degrees. Although, the projective data association can be imitated with LiDAR data by projecting the point cloud within the local sliding window onto a spherical image plane, this causes uneven spatial resolution where the top and bottom is dense and around the equator is sparse. Dividing the spherical image plane into an equal gird for an even sensing plane requires a nearest neighbour search, significantly reducing the efficiency of the projective data association method.
3) Map Element Fusion: Once matching of map elements are established, the next step is to fuse the measurements to improve estimation. Keller et al. [24] proposed a dense surfel fusion method by simplifying the Bayesian estimation from 3D to 1D. In their approach, each map element is independently updated, making its computation much simpler than the EKF case. They also utilized radial distortion as an initial uncertainty parameter and reduced the uncertainty whenever the surfel is observed again. ElasticFusion [1] further extended the uncertainty as a function of sensor motion to consider the uncertainty caused by motion blur. However, those simplified Bayesian models are not appropriate for dense 3D LiDAR mapping where the existence of a surfel degeneracy often causes slower convergence.
The approaches in this category mostly do not consider the sensor noise model [1], [6], [23], [24] or they are dedicated for RGB-D models only [26].

III. OVERVIEW
Our proposed map-centric SLAM system is composed of three main components: local mapping, global mapping and loop closure detection as shown in Figure 3. The local mapping part takes visual, IMU and LiDAR measurements to build motion-distortion-corrected maps using continuoustime trajectory optimization. This stage is similar to the sliding windows in [3], [21], [33], but it is different in that it takes the sparse global map as a map prior. This let the new local map always registered to the global map.The loop closure detection module keeps generating 2D features from the visual sensor and compares them to the previously generated key frames.
The global mapping part builds and fuses surfel maps. Note that two different surfel maps, multi-resolution 3D sparse surfel map and dense surfel map, are utilized for different Continuous-time trajectory representation The indexes k of local control points to recover   [21] are ideal for fast and robust continuous-time trajectory optimization, whereas it is too sparse for dense representation [22]. Thus, we utilize a fixed size flat-shape hexagonal surfels [24], [27], [34] for dense surfel fusion. Global map consistency is achieved by non-rigid deformation using the estimated misalignment from the loop closure module. It is inspired by the global deformation in ElasticFusion [34], and we further extend it by considering surfel uncertainty propagation in LiDAR. Upon a visual loop closure detection, the loop closure module finds the 6 DoF misalignment between the global map and the current frame over different places until the uncertainty of the misalignment is below a certain threshold. Table I compiles all states and parameters used in this work.
There are two types of surfel map representations used in this paper, sparse surfel map S g and dense surfel map M g . Each surfel map is individually updated with their local maps S l and M l which are extracted from the current laser scan.
The sparse surfel map consists of 3D ellipsoids extracted from laser points using multi-resolution voxel hashing [21]. Each ellipsoid is defined with a centroid c ∈ R 3 and a covariance matrix Σ c ∈ R 3×3 which represent the distribution of points within the voxel. Likewise, the dense surfel map maintains the mean and covariance defining the ellipsoid of points in the surfel, which in turn allow fused estimates of position p ∈ R 3 and normal vectorˆ n ∈ R 3 . In contrast to conventional surfels [1], [24], we associate uncertainty Σ p ∈ R 3×3 and scatter Ξ p ∈ R 3×3 with the position of each disc surfel, which are later used to merge surfels based on Bayesian filtering. The normaln of surfels are extracted from the scatter matrix Ξ p . Note that 3D ellipsoid surfels are expressed with ellipsoids of their covariance matrices, while 2D disc surfels are expressed with discs with normal directions.

IV. LOCAL MAPPING
To deal with the asynchronous estimation and LiDAR motion distortion, we are introducing the concept of the continuous-time trajectory representation and optimization to the map-centric approach. In this section, we describe our proposed constraints for the map-centric local trajectory optimization that systematically couples the local continuous-time trajectory optimization method with the map-centric method. Also, in the later part of this section, we describe our strategy that improves the accuracy and efficiency of the trajectory optimization and update.

A. Continuous-Time Trajectory Representation
Let T(τ ) be composed of translational component t(τ ) ∈ R 3 and rotational component R(τ ) ∈ SO(3) as, Then, utilizing the linear continuous-time trajectory representation, its value can be evaluated by an interpolation from two poses T k , T k+1 where their timestamps satisfy τ k < τ < τ k+1 . Given poses with timestamps, their interpolation at τ is given by, Fig. 4: Illustration of geometrical constraints on the local trajectory. Two surfels u τa and u τ b generated at τ a , τ b within a local window forms constraints on the interpolated trajectory poses T τa and T τ b . The constraint between the local scan u τc and the map prior u m forces the trajectory to be fitted onto the map prior.  where the relative pose [ξ] × ∈ se(3) is defined by log(T −1 k T k+1 ) and the exponential mapping e α[ξ]× linearly interpolates the relative pose on the manifold with the inter-

Map Prior
The linear interpolation error proportionally increases when the sensor is in high-speed motion. The error is propagated to the map deteriorating the map quality. To reduce the effect of the linear interpolation error, we have generated the trajectory at 100Hz.

B. Local Trajectory Constraints
When the motion of a device is relatively moderate compared to the scanning speed and the initial trajectory guess is fairly good, shapes of local features are often well preserved in the measurements even with the motion distortion. The first two geometrical constraints utilize this property of LiDAR scans [21], [35]. This process starts by transforming the LiDAR measurements with respect to the world frame and extracting sparse ellipsoidal surfels from multi-resolution voxels [21]. Then, we calculate the point-to-plane errors between two corresponding new surfels a and b in their averaged normal directions n ab as where u τ is the centroid of a surfel, t(τ ), R(τ ) are the interpolated sensor pose at time τ from the continuous-time trajectory. The visualization of the constraint is given in Figure  4 and Figure 5.
The second geometrical constraint defines the relationship between the current local map and the map prior. The map prior is obtained from the active global map which we will be described in the Session V. The map prior constraint makes the local trajectory always aligned to the previously built global map, which is essential for map-centric operation. Thus, we define the constraint between a new surfel c and its corresponding global map surfel m as where n mc is the averaged normal directions of surfel m and c. Note that we do not have to transform the map prior u m because they are points in the world coordinate system. The motion distortion is corrected based on the map prior. The initial map prior is given from a short period of stationary scanning at commencement. On the other hand, the inertial information from IMU offers a prior on the trajectory in terms of rotational velocity and translational acceleration. The IMU measurement constraints are written as, where α, are respectively acceleration and rotational velocity from IMU measured at time τ , g is the gravitational acceleration, R(τ ), t(τ ), ω(τ ) are respectively interpolated rotation, translation, and rotational velocity of the continuoustime trajectory at time τ , and b α , b ω are the IMU biases. These two constraints form a strong restriction on local smoothness.

C. Trajectory Optimization and Update
Compared to the number of the interpolations required for the motion distortion correction of the point cloud, the number of the interpolation for the trajectory update is much fewer. It is not necessary to compromise accuracy for a fast interpolation by using linear interpolation for this stage. Therefore, we utilize B-Spline for the smooth update.
We start by defining a set of control points Q k ∈ Q which is created at an interval of τ step and longer than the discrete trajectory sampling interval. The trajectory corrections are achieved by finding a set of parameters that minimize the residuals of Eq. (3)-(6) as, where . d is a time lag parameter and Q is control points. Then, the correction to the discrete trajectory is made by, where the correction dT for each discrete trajectory at τ k is defined as an interpolation of neighboring trajectory elements represents a set of neighbor indices at τ . Although the utilized control points are in Euclidean space as c rk , c tk ∈ R 3 , in each iteration they always starts from zero. Thus, it is free from the singularity problem.
For the recovery of the correction dT from the control points we have utilized the B-spline as, where t ∈ [0, 1]. In a similar way, rotational component r τ ∈ so(3) is defined from separate control points c r and converted to SO(3) by an exponential mapping R τ = e [rτ ]× . The approximation method in [3] utilized SO(3) + R 3 update scheme where the translation update is given as (t = δt + t). The selection of SO(3) + R 3 over SE(3) was reasonable as the SE(3) update with (t = δt + e [ω]× t) could cause an accuracy problem with linearization. Especially when t is large, the linearization error of the rotation enlarges the translation error in the second term e [ω]× t. Separately updating SO(3) + R 3 reduces the linearization error in the large map global optimization. However, SO(3) + R 3 results in suboptimal. Considering that our approach does not need a global trajectory optimization, we use a SE(3) update. Also, as the map grows, this problem can be solved by simply converting the trajectory and estimations in the local window to be relative to the first frame in the local window.
After the optimization, local dense and sparse surfel maps M l , S l are built and updated by the optimized trajectory. Each dense surfel in the local dense map is composed of position p ∈ R 3 , surfel normal n ∈ R 3 , and timestamp t. Also, the uncertainties of position and scatter matrix Σ p ∈ R 3×3 , Ξ ∈ R 3×3 , which are utilized in dense surfel matching and fusion are calculated from their neighboring points. On the other hand, local sparse surfel maps S l which are utilized in Equation (3), (4) are updated by the optimized trajectory. Each sparse surfel is defined with a centroid c ∈ R 3 and a covariance matrix Σ c ∈ R 3×3 which represents the distribution of points within the voxel. The local maps M l , S l are fused into the global maps M g , S g in the following section.

V. DENSE SURFEL MATCHING AND FUSION
In this section, we introduce our probabilistic surfel fusion. We will detail our two staged probabilistic surfel matching that can maintain surface resolution of the dense map to a specified value. Also, a new probabilistic surfel representation that can correctly track the surfel normal under the existence of degeneracy will be highlighted.
A. Surfel Uncertainty Modeling 1) Uncertainty in Position and Orientation: We model the surfel position and shape as a random ellipse using a normal inverse Wishart model, based on the statistical model of LiDAR sensor noise measurements [36]. If the LiDAR points were noise-free, the normal inverse Wishart model would provide an exact closed-form estimate of the mean and covariance matrix from which the points were drawn. In practice, the positional uncertainty of the LiDAR points is higher along the beam direction. Therefore, we adopt the method of [37], which provides estimates under nonhomogeneous noise matrices. By estimating both the mean and covariance matrix, we recursively integrate knowledge on both the position and orientation of the surfel.
Surfels are modelled as a Gaussian distribution of points with mean and covariance to be estimated. The LiDAR uncertainty is an additional Gaussian measurement noise, such that the distribution of each point is N {z; µ, Σ + Q w }, where µ and Σ are the unknown surfel mean and covariance, and Q w is the measurement noise in world coordinates. Since points in the same surfel that are local in time will have the same geometry, this can be modelled as where w R l and l R b are rotation matrices from laser to world coordinates, and from beam to laser coordinates, respectively. The amount of uncertainty along the beam direction is defined by the sum of the distance uncertainty σ 2 r and additional uncertainty σ 2 i caused by the incident angle [38], [39]. The uncertainty covariance in beam coordinates is where σ 2 r is noise due to the beam radius, σ 2 d is the nominal depth variance following the model of [36], and σ 2 i is an additional uncertainty caused by the incidence angle [38], [39]. Surfel observations from different times are assumed conditionally independent given the surfel mean and covariance.
2) Uncertainty in Colour: The local point cloud is rendered to the image frame for visualization. We parameterize three uncertainty factors that affect the accuracy of the rendered point cloud as α r = e(r/r th ) − 0.5, attributes in a sigmoid function, the colour uncertainty σ c is given by,

B. Surfel Matching
The local surfels M l are transformed into world coordinates to find matched surfels in the global dense surfel map M g . The matching process begins with finding a set of candidate surfels A g for each surfel ϕ l ∈ M l . For efficient matching, initial matching candidates are selected by the octree-based nearest neighbor search algorithm.
Then, the resolutional distance r between each source and destination pair in Figure 6 is compared with a resolution threshold θ r to decide if their projections are close enough on the surface. If so, we check the depth d in the Mahalanobis distance. To consider the uncertainty only along the normal direction, we propagate the positional 3D uncertainties of source and destination surfel Σ d , Σ s into 1D along each normal direction by σ 2 =n Σn.
Finally, if the 1D Mahalanobis distance along the surface normal direction is less than a threshold θ d , we assume that they are in correspondence, and put the matched surfel into B g . Note that the resolutional distance d is compared in the Euclidean space to preserve the desired surface resolution in the Euclidean space. Algorithm 1 summarizes this surfel matching process. Note that our matching method enables the matching process to search more along the beam direction, while effectively maintaining the desired surface resolution without a voxel grid.
The RGB values of each matched point in the image are found by rendering the surfel map at each camera position.

C. Surfel Fusion
As described previously, surfels are modelled as a Gaussian distribution of points with mean µ and covariance Σ. These

Algorithm 1: Finding Surfel Matches
Input: Global dense surfel map M g and a surfel ϕ l ∈ M l Output: end 10 end parameters are estimated recursively using a normal inverse Wishart model. In each update, the LiDAR points {z i } are summarised via the noise uncertainty Q w given in (10), alongside the number of points n, the measurement mean and the scatter matrix:z As described further in Appendix A, the method in Algorithm 2 recursively updates the surfel parameters, based on the random matrix approach in [37]. By directly estimating the surfel ellipsoids, surface centroids and normals are estimated through an integrated Bayesian model. Normal estimates may then be extracted from the estimated ellipsoid covariance via eigen decomposition. This is different from our previous work [2], [10] and other approaches [26], where normal vectors were extracted from each surfel, and fused by assuming an error Normal Extraction: Colour Fusion:

end
model on these estimates. Our new approach uses an error model that is implied by the single point error model. The colours are fused similarly with a Gaussian distribution but in R as there is no measurement to estimate the uncertainty of each colour channel. The surfels that are not matched to the global map will be added to the global map as a new unstable surfel. To effectively remove the surfels generated from LiDAR data with non-Gaussian noise, any unstable surfels that are not re-observed on revisits will be deleted from the global map.

D. Active and Inactive Maps
One of the most important assumptions in the fusion is that the global map consistency should be guaranteed around the area where a fusion occurs. In the worst case, surfels will be fused right before the loop closure. The fusion of the surfels with a misalignment will make an unwanted deformation in the local area. Thus, similar to [34], we introduce active and inactive maps A, I according to the surfel timestamps so that new surfels are always optimized and fused within the active area. To detect the misalignment between the active and inactive areas, it finds the amount of overlap and misalignment by rigid Iterative Closest Point (ICP). For a robust and fast ICP, only sparse surfels are utilized with a geometrical weight. When the overlap and misalignment is large enough after the ICP, it triggers deformation that will be described in the next section. To maintain the map coherency between active and inactive components, matched inactive components to active components are found and updated in every step. Algorithm 3 shows our temporal fusion method for LiDAR surfels.

VI. GLOBAL MAP BUILDING
As the proposed system does not maintain any trajectory, an alternative method for a loop closure is necessary. In this section, we describe our approach for maintaining the map to be globally consistent where we elastically deform the map upon the loop closure detection. Also, we highlight our sequential metric localization method that can stably and robustly estimate the misalignment.

Algorithm 3: Temporal Surfel Fusion
Input: New Frame S l , M l , Global Maps S g , M g Output:

A. Map Deformation
Upon loop closure detection, map deformation is carried out to maintain global map consistency. Deformation nodes and loop closure constraints are only selected from the sparse surfel map S g . Then, once the optimal deformation is found by an optimization, it is applied to both entire map elements in M g , S g . We adopted the graph deformation technique for the introduction of elasticity in the map and temporally connected deformation graph of [34]. However, our approach is different from [34] in that the number of deformation nodes is decided by the area of the reconstructed space and the uncertainty is also deformed along with the normal and centroid.
1) Graph Construction: Graph nodes are constructed from randomly selected surfel centroids of S g . The centroids are selected to uniformly represent the space. The number of nodes is decided by the area of the space. As the proposed dense surfel map fusion maintains a canonical form of surface without redundancy, the area of the map surface, A can be easily found by counting the number of total surfels. Thus, considering that each surfel of the map has a fixed radius circle, the total number of nodes n is, where γ is the total number of map elements, r is the radius of the circle, and E is the number of nodes in the area. Denser nodes help reduce local spatial distortion but exponentially increase the computational cost. To reject any temporally uncorrelated connection, we order the nodes by temporal sequence and define their neighbors according to their temporal closeness. The surfel timestamps are updated to the current time when they are revisited. The deformation graph is composed of node positions g j ∈ R 3 , node rotations and translations, R j , t j ∈ SE(3), and a set of neighbors V(g j ). Most of the time, the map is deformed by the difference between node translations. Regarding the number of neighbors, previous research [34] empirically indicates that four neighbor nodes are sufficient.
2) Graph Deformation: Assuming that a set of graph node locations g j are established and their deformation attributes R j , t j are found, the graph deformation can be applied to the entire map. Given a set of node parameters, an influence function φ(p i ) that deforms any given point p i by the j-th node g j is defined as, For a smooth blending of deformation, the complete deformation of a point is defined as a weighted sum of the influence function φ(p i ) with its neighbor nodes N(p i ) around the position. When selecting nodes, their temporal and spatial closeness are considered. Thus, the deformed pose p i is defined as, where the weight w j (p i ) is decided by the distance between the node g j and the point p i , where d max represents the max distance between the node and point within N(p i ). Those nodes that are relatively far from the given point p i have smaller weights and less effect. Note that Equation (16) will be applied to both sparse and dense maps S g , M g . Not only surfel centroids, but other surfel attributes also need to be deformed. The new normal direction of the surfel is defined as, The uncertainty and geometry matrices Σ pi , Ξ i , Σ ci are deformed by first order linear propagation. The general form of the propagated uncertainty and geometry in the new deformed space is given as, where R i represents the blended rotation of the surfel location and defined as, 3) Graph Optimization: In this section, we will describe constraints for the graph optimization, given a set of matched surfel centroid pairs p src , p dest ∈ S g . During the optimization, the transformation R j , t j of each node which minimize the source and destination will be found. The first constraint that reduces the difference between the deformed source sets p src = φ(p src ) and target surfel p dest is, The distortion is applied over all the global map surfels including the destination surfels p dest itself. To prevent the infinite loop of deformation between source and destination, the following pinning constraint will fix the destination surfels not to be deformed, To guarantee the smooth deformation over the whole region, the following term spreads the deformation to the neighboring nodes. When the node rotation R j is near the identity matrix, the smoothness term indicates the distance between t j and its neighboring node translation t k , (23) Then, the optimal node transformationR j ,t j that minimizes radical deformation and loop closing error is defined as, where ω reg , ω pin , ω loop represent the weights for each constraint to tune elastic property. We follow the proposed values in [34].
Since the cost function is a non-linear pose graph optimization problem in the form of f (T) = Tp + c on a manifold, we use the non-linear iterative Gauss-Newton method.
4) Global Loop Closure: There are two sources of loop closure detection. For a loop where the source and destination distance is moderately close, the ICP in Algorithm 3 will detect the misalignment first. However, for large-scale mapping where large misalignment can occur, an alternative place recognition such as [40], [41] is required. We propose a robust method for place recognition and metric localization using both visual and geometrical information, which will be detailed in the following section. In both cases, loop closure constraints are given as follows, where p src are randomly selected points from S l and R, t denote misalignment between S l and I. Note that we utilize a transformed destination p dest from the p src to prevent any unwanted deformation caused by source and destination point difference. This way, the deformation implicitly reduces the misalignment [R, t] to be [I, 0].

5) Misalignment Estimation for Global Loop Closure:
Disadvantage of the map-centric method is that the result of an incorrect loop closure is more destructive and is not reversible. Thus, for reliable estimation of the misalignment [R, t] even with incorrect initial guesses and large drifts, we tightly fuse the surfel-based point-to-plane constraint with the 3D feature based point-to-point constraint as, where the pair of the matched surfel centroids p r , p s are acquired by the nearest neighbour method in the surface normal and centroid space. The matched sparse 3D features p r , p s are found through 3D feature matching such as FPFH [42], but not limited by the types of feature descriptor. The TABLE II: Comparison of different continuous-time trajectory optimization in a simulation. Our method has significantly lower state dimension, compared to the same motion resolution in [9] but also higher accuracy than [3].
surfel based constraint and the 3D feature based constraint are complementary to each other. While the surfel constraint lacks semantic correspondence which result in sub-optimal solution, the 3D feature often suffers from noisy LiDAR point cloud or partial observation. When combined, the 3D features roughly guide the optimization to where the surfel correspondence can be correctly found and surfel constraints refine the final result. Similar to our previous approach [11], [43], we sequentially estimate [R, t] at different places over time which greatly reduces the chance of failure. Refers to the Appendix B for the details of the sequential pose fusion.

VII. EXPERIMENTS
In this section, we present the quantitative and qualitative performance analysis of the proposed method in terms of trajectory estimation and reconstruction quality with various sensor payload configurations and environments. Also, we provide a qualitative comparison of the metric localization accuracy with various approaches.
We have utilized CT-SLAM [3] as a baseline method, which was first published in 2012, but has been optimized over the last ten years proving its robustness and accuracy under various robotic platforms and environments [44], [45].

A. Simulation on Trajectory Optimization
We demonstrate the performance of the proposed method by comparing it to the most popular spline based trajectory optimization model through simulation in order to have quantitative analysis with the trajectory ground truth.
We have utilized a five second length of local window with a simulated angular velocity and linear acceleration at 100 Hz. The initial trajectory t 0 is built by an accumulation of the simulated angular velocity and linear acceleration with a bias and Gaussian noise. Randomly generated sparse surfel features (1000) are utilized for the trajectory optimization and assigned with a timestamp. For simplicity, we have only utilized the local surfel to fix the surfel registration constraint in Eq (4). The simulated trajectory optimization and corrections are visualized in Figure 7.
To demonstrate the advantage of the proposed method in local trajectory estimation, we present the comparison of the composition trajectory model and approximation model along with a few different parameters in Table II. Interpolation represents the interpolation method utilized for the cost function. The result in Table II shows that when the update is smoothed by spline and properly done on SE(3) the Continuous-Time trajectory optimization result is better than the linear SO(3)+R 3 composition method originally proposed in [3] and the approximation model [9].

B. Experimental Setup
To demonstrate utility of the proposed method under various scenarios, we have collected multiple datasets with different sensor setups mounted to multiple robotics platform as in Figure 8. The collected datasets are categorized into three  classes according to sensor types, environments, and robotic platforms.
Through Figure 9 (a) to (c), the datasets are collected with a single beam hand-held 3D spinning LiDAR. The device consists of a spinning Hokuyo UTM-30LX laser, an encoder, a Microstrain 3DM-GX3 IMU, RGB camera. In the collected datasets sensors were moved with 0.9 m/s and 0.7 rad/s while rotor spins at 1 rotation/s. The datasets in Figure 9 (d) to (e) are collected with a hand-held multi-beam LiDAR with Velodyne VLP-16, a Microstrain 3DM-GX3 IM and an independent Gopro camera. Gopro does not share a common clock with LiDAR.
The dataset in Figure 10 (a), (b) are presented to demonstrate the utility of the proposed method on realistic robotic platforms, respectively, a four legged robot and ground robot. For those dataset, the same sensor payload (Velodyne VLP-16 LiDAR, IMU and camera) was utilized while they are mounted on each platform as in Figure 8.
Also, we have collected datasets in different environments to demonstrate the robustness of the method against the different geometrical shapes. The datasets (a) and (b) in Figure 9 are collected indoors whereas Figure 9 (c) to (e) are collected outdoors. The outdoor dataset mixes structured (c) and (d) and unstructured places (e).  Spatiotemporal extrinsic parameters between the LiDAR and the rest of the sensors (camera, IMU) are estimated using the method in our previous publication [46]- [48]. The spatiotemporal parameters are utilized for a colorization of surfel. Shader-based surfel render is developed for the colorization and 3D scene rendering. Some of the scenes of the dataset are rendered and presented in the third row of Figure  9. For the colorization, a distorted image is directly utilized along with the calibration parameters in the shader to avoid loosing information and reduce computational cost. Some of the datasets require a global loop closure detection. The global loop closure is found by 3D and 2D combined loop closure detection method [11], [45] and the misalignment is estimated using the sequential method described in Section VI using only 3D point cloud. Detected loops are closed by the deformation graph as depicted in Figure 11 (b).

C. Trajectory and Structure Comparison
To evaluate the trajectory estimation accuracy of the proposed method, we compared the deformed trajectory with the globally optimized trajectory from the benchmark method [3]. Note that in the proposed method, the trajectory is not required to be stored. They are saved and deformed only to evaluate accuracy of the proposed method. For the comparison, we utilized the absolute trajectory Root Mean Square Error (RMSE) metric, which estimates the Euclidean distances between the deformed trajectory and the globally optimized trajectory. Table III shows the error and trajectory statistics of each dataset visualized in Figure 9. The trajectory error of each estimations varies from 0.2 to 0.5 meters with high correlation to the map size. The trends shows that the difference gets larger when the map scale is growing, which is reasonable considering that the proposed method does not count gravity direction in the deformation graph optimization. The initial position is aligned with the gravity direction but when the sensor scans further area from the initial location the gravity alignment faints, making the difference. A qualitative comparison of trajectories from the dataset (b) is presented in Figure 12. The visual comparison shows that local details are well preserved and similar to the globally optimized trajectory, but there were overall 0.1 meters offset between the two trajectories.
Also, to visualize the structural difference of highly redundant scanning case, we registered the two point clouds respectively from the baseline method [3] and the proposed method and calculated the point-to-mesh distances in Figure  13. For the point-to-mesh distance, a fine mesh is extracted from the globally optimized point cloud of the baseline method [3]. The result shows the reconstruction difference was within ±0.02 meters but it was growing on the outer side.

D. Surface Estimation
For a quantitative comparison on real datasets where ground truth is not available, known planar surfaces such as a floors or walls are utilized for calculating the relative noise level. Table IV shows statistical analysis of the multiple planar areas. The error in the table represents the mean projective distance, which is calculated from the mean plane of each patch. With the same data extracted from the planar patches, the analysis result on the relationship between the number of fusion steps and accuracy is given in Figure 14 along with the comparison with the original raw points cloud on the left most side of each figure. For relative comparison, the distance of each points in the area from the average plane of the area were used to calculate normal direction errors and projective distance errors. γ represents the number of fusion where surfels with the same number of fusions are utilized for the statistics. For example, for the calculation of γ = 2 surfels that are fused only 2 times are utilized. The result in Figure 14 indicates that the errors drastically are reduced when surfels are fused more than four times reaching ±10mm   and 0.2rad (≈5deg) error range. The surfels with low numbers of observations are removed from the global map after some period of time as they tend to have higher possibility of being outlier or mixed pixels [49]. The development of the map and the surfel statistics changes over time are visualized in Figure  15 with the dataset (f). The figure shows how the number of surfel observations impact uncertainties when the sensor makes multiple traverses on the same place.

E. Evaluation of Localization for Loop Closure
For the initial loop closure detection, we have utilized the visual place voting method [40]. Once a place recognition triggers the proposed localization for the loop closure, a procedure for estimating 6 DoF misalignment starts where it finds only single misalignment between the very first revisited place found from the previous map and the new map. In the sense that the process includes 6 DoF registration with unknown initial guess, the problem is closer to a large scale global registration problem [50] or kidnap recovery problem rather than a place recognition problem [51]. To this end, we compare our method with the well known global ICP registration libraries [52], [53] with different initial guesses.
We have utilized outdoor and indoor mixed point cloud datasets to evaluate the localization performance. The globally optimized trajectory was utilized for calculating the ground truth. To estimate the robustness against initial guess, initial poses are randomly generated with ten different loop closure scenarios. We have set the covariance for the random misalignment generation to be larger along the z axis (gravity direction) and smaller in x and y axis (orthogonal to gravity direction). The 6 DoF misalignment is estimated 50 times each with ten different locations which simulates 500 loop closure triggers. For each loop closure trigger, misalignment starts sequential estimation over different places along the trajectory which exactly simulate the loop closure procedure with real data.
The RMSE of the estimated pose is listed in Table V along with sparse ICP algorithm (a), and the proposed method (e) and the other global registration algorithms (b), (c). In (b), Open3D [52] global registration is utilized with an initialization by FPFH [42] feature and RANSAC. Similar to (b), SHOT was utilized for the initial pose estimation, then refined by point-to-plane ICP in (c). Longer localization processing time drastically increases map fusion cost between the active map and the global map. Hence, we are not comparing to the time consuming algorithms such as GO-ICP [50] which could take several minutes to hours to process a single pair of point cloud.
The experiment shows that the surfel ICP (a) is prone to a local minima problem with a large displacement. The  large displacement frequently occurs especially in a large scale mapping scenarios. The pose estimation from Open3D registrations (b) is generally reasonable regardless of the initial guess. However, occasional registration failures occurred due to the failure in finding feature matching. SHOT (c) showed a moderate rotation estimation, yet the estimated translation often suffered from outliers even with a close initial guess. This is assumed to be due to the failure in the feature based pose initialization. The proposed method (d) shows the most stable result. Although, both translation and rotation estimation is slightly increasing reaching up to 6 cm and 0.004 rad (≈0.22 deg) in the hard case scenario but they are within a reasonable range and still better than the other algorithms.

A. Trajectory Optimization
The composition model stores 500 raw poses and only optimizes the trajectory corrections whereas the approximation model entirely represents the trajectory as the abstracted few spline control points (Figure 7). Because of this difference, with the same state dimension to be optimized, the approximation method has significantly lower accuracy. With relatively higher motion resolution (101 control points), the approximation method follows the ground truth trajectory reasonably but high frequency motions are ignored. The ignorance of the high frequency motion propagates noise to the map when the scanning equipment moves in a high frequency motion such as in the case where the device is mounted on a wheeled platform driving outdoor. Considering LiDAR measurements reach up to several hundred meters, even a small angular difference can cause a large displacement of the points.

B. Map Estimation
While the local details are very similar to the baseline method in Figure 12, there were overall structure movement about 0.1 meter which make similar movement in the deformed trajectory as shown in Figure 12 (b). This could be problematic when the scanning scenario includes partial revisit CT-SLAM [3] Proposed method   For each noise level, initial poses are randomly generated according to the following parameters: Easy σ θz=10 , σ θxy=1 , σ t=0.5 , Medium σ θz=50 , σ θxy=5 , σ t=5 , Hard σ θz=100 , σ θxy=20 , σ t=50 (σ θ in deg, σ t in meter).
of a space which does not make enough overlap to trigger a local loop closure. While the map-centric approach solves the time complexity problem, this is fundamental weakness of the proposed method when used with a long range LiDAR. Thus, we suggest to reduce the fusion range.
The highly redundant scanning case in Figure 13 showed slightly different pattern to the dataset (b). The center part is close to the baseline method but outer part of the map showed difference of ± 0.02 meters. This is due to the loop closure occur in a the first traverse as visualized at (i) in Figure 11 (a). The short distance make relatively small misalignment at the loop closure therefore the deformed map around this part is similar to the baseline method. However, when a new area is explored and the map grows, the difference starts to increase. Note that the difference in this case is not because of the map deformation but because of the accumulated small drifts.
Statistics in Table IV shows that the estimated surface is up to three times less noisy than the point cloud generated using baseline CT-SLAM [3]. However, repeated non-Gaussian noises such as noises from mixed pixel problems [49] are not effectively fused by the proposed method. The repeated non-Gaussian noise pattern usually appears from stationary scans.
Regardless of the number of the original raw points in the proposed method, the surfel density in patches are uniformly maintained after fusion. Figure 16 visualizes a small patch of the dataset (e) with a reduced surfel size to show that the surface forming surfels are non-redundant and uniformly distributed. However, because of the restriction in surface resolution, objects which are smaller than the given surface resolution are often ignored, losing some details of scene. When, an excessive surface resolution is utilized computational cost significantly increases due to the matching and fusion cost and also without increasing the reconstruction quality. With our experimental sensing payloads, we found that 0.02 meters or higher resolution were acceptable. The result in Figure 15 shows that after the first traverse of the map (t=100 sec), the number of observations of each surfel increases (second column) while the uncertainties decrease (third and forth column). The uncertainty of the overall surfel positions decrease whereas the scatterness along the z direction of the surfels around edges does not decrease after some point. This is due to the way the scatter matrix is extracted from the point cloud where surfels around edges and corners always have a high eigenvalue in the z direction.

C. Localization for Loop Closure
Compared to the methods (a), (b), (c) in Table V our proposed sequential metric localization method utilizes more point clouds collected at different locations. However, rather than utilizing a single large chunk of the point cloud, the proposed method (d) segments the map into frame level and runs independent registrations. While this approach provides an evaluation metric to decide the reliability of the current estimation, it also offers cross-validation ability between the independent misalignment estimations. These two components in the proposed method significantly contributed in reducing localization failures. However, the translation error is still relatively higher than the map reconstruction accuracy. Considering the surface reconstruction quality is often goes below 1 centimeter, it is desirable to reduce the metric localization accuracy to the similar level.

IX. CONCLUSION
We presented a new approach for dense LiDAR-based mapcentric SLAM. Our main contribution in this paper is the identification of the problems that occur when applying mapcentric dense mapping approach to LiDAR sensors, followed by a new SLAM framework that successfully adapted the beneficial characteristics of the map-centric approach and continuous-time representation.
We have tested our method with 8 dataset collected under various challenging places in different scales with 4 different sensor types and configurations. Our experimental results indicate that, the deformed map showed only 2 centimeters difference, while the surface estimation is up to 3 times less noisy after multiple revisits. Also, the sequential misalignment method gives a stable and accurate localization estimation within 6 centimeters even with the hardest cases where other common methods fail.
While we saw the possibility of LiDAR device for mapcentric SLAM, we found couple of limitations such as growing map distortion or the partial observation problem. The distortion in the 60 meter scale map reached to 10 centimeter and expected to be larger with extended map size. However, we believe these weakness can be solved by further introducing additional priors and constraints. So, we leave these problems as future works.

APPENDIX A SURFEL FUSION
Here we briefly describe the method that we utilise for fusion of surfel ellipsoids. If points were observed noise-free, i.e., LiDAR points were distributed according to N {µ, X} where µ and X are the true surfel parameters, then a normal inverse Wishart model may be used to estimate the implied ellipsoid by accumulating the overall mean and sample covariance. A fast, approximate method for accommodating noisy observations of the form N {µ, X + Q} is proposed in [37]. The method maintains a state which includes the estimate of the surfel centroidμ, the covariance of this estimate,Σ, the accrued sample matrix for the extent,Ξ, and the count of points included in the estimates, ν. The updates due to a new surfel with n points with meanz, scatter matrixZ and measurement noise Q w are: S =Σ + Y n Y =X + Q w (32) where n z is the measurement dimension, an the extent update of: N =X 1/2 S −1/2 NS −T /2X T /2 (36) The estimate of the covariance defining the object ellipsoid iŝ X. An alternative method using variational inference has also been proposed [54], but this requires iterative solution, and is avoided for performance reasons.

APPENDIX B SEQUENTIAL POSE FUSION
The Bayesian fusion provides a closed-form solution on the vector fusion problem [2]. However, directly applying the Euclidean Bayesian fusion on the poses causes a convergence to a suboptimal as the pose vector is on manifold. Thus, we utilize the sequential SE(3) pose fusion approach proposed in [11], [55].
This begins by modeling an error n ∼ N (0, Σ n ) between an individual estimation of a pose T n and the ground truth where the second and third equation substitute the ground truth which is not directly available and convert the error as a function of a small perturbation ξ ∈ se(3) regarding our best guess so far on the pose T * and the individual estimation. Thus, given N poses the best pose that minimize the sum of the error can be achieved by iteratively minimizing the following cost function regarding ξ.
where I −1 in Equation (41) is Baker-Campbell-Hausdorff approximation [55] of Equation (39) which is also referred to the left jacobian of SE(3). The pose fusion model above suggests a framework for the SE(3) pose fusion problem. The multiple misalignment estimations in different places are fused by this framework. However instead of following the original batch fusion, we sequentially fuse posed until the uncertainty is small enough.
Let T c , Σ c , T k , Σ k , respectively be k th new alignment estimation and fused alignment measurement up to the point with their corresponding uncertainties. Note that ideally T c and T k are identical as the new alignment estimation T k is always represented with respect to the first frame where the first alignment was estimated.