{'fraction_non_alphanumeric': 0.03127170524479478, 'fraction_numerical': 0.011145390815551858, 'mean_word_length': 4.616947922337144, 'pattern_counts': {'":': 0, '<': 0, '<?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': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
LiDAR Road-Atlas: An Efficient Map Representation for General 3D Urban Environment

In this work, we propose the LiDAR Road-Atlas, a compactable and efficient 3D map representation, for autonomous robot or vehicle navigation in general urban environment. The LiDAR Road-Atlas can be generated by an online mapping framework based on incrementally merging local 2D occupancy grid maps (2D-OGM). Specifically, the contributions of our LiDAR Road-Atlas representation are threefold. First, we solve the challenging problem of creating local 2D-OGM in non-structured urban scenes based on a real-time delimitation of traversable and curb regions in LiDAR point cloud. Second, we achieve accurate 3D mapping in multiple-layer urban road scenarios by a probabilistic fusion scheme. Third, we achieve very efficient 3D map representation of general environment thanks to the automatic local-OGM induced traversable-region labeling and a sparse probabilistic local point-cloud encoding. Given the LiDAR Road-Atlas, one can achieve accurate vehicle localization, path planning and some other tasks. Our map representation is insensitive to dynamic objects which can be filtered out in the resulting map based on a probabilistic fusion. Empirically, we compare our map representation with a couple of popular map representation methods in robotics and autonomous driving societies, and our map representation is more favorable in terms of efficiency, scalability and compactness. In addition, we also evaluate localization accuracy extensively given the created LiDAR Road-Atlas representations on several public benchmark datasets. With a 16-channel LiDAR sensor, our method achieves an average global localization errors of 0.26m (translation) and 1.07 degrees (rotation) on the Apollo dataset, and 0.89m (translation) and 1.29 degrees (rotation) on the MulRan dataset, respectively, at 10Hz, which validates the promising performance of our map representation for autonomous driving.


Introduction
Now days, autonomous online 3D mapping plays a key role on robot exploration of unknown environment, and can also acts as the "teaching" role in most "teach-and-repeat" paradigms for robot navigation or autonomous driving applications. As exploration goes on, the map size (either runtime memory usage or hard-drive storage) can scale up quickly. Likewise, 3D mapping of large general urban environment is also important for commercialized autonomous vehicle, where the resulting 3D high-definition (HD) map can be re-used for vehicle localization and path planning. When applied at large scale, the map's storage size is one key factor that has to be considered as well. Therefore, a parsimonious map representation is highly desirable. In this paper, we propose the LiDAR Road-Atlas, a compactable and efficient map representation of general 3D urban environment, for autonomous robot or vehicle navigation. It will be generated by an online mapping framework based on incrementally merging 2D local occupancy grid maps.
In the robotics society, one classical and efficient way to represent the environment is the 2D occupancy grid map (2D-OGM) [Elfes, 1987] where the cells in the grid are labeled as drivable, unknown or occupied. A 2D occupancy map usually suffices for robot navigation in structured outdoor environment. For navigation in non-structured outdoor environments, where ground is not purely flat ( Fig.1) or where overhanging objects (e.g. tree branches, tunnels, low roofs, etc.) are very common, a 2D representation of the environment is insufficient because 2D-OGM cannot fully represent the actual environment. In such outdoor scenarios, 2D-OGM cannot be accurately created in general because it is not easy to reliably find the boundaries between known (traversable region) and unknown (occluded region) areas due to the existence of low obstacle regions. Using an underrepresented map can be totally erroneous in path planning, let alone vehicle localization. Therefore, a more suitable alternative for outdoor navigation would be a map representation where the environment could be well described. Over the years, a few more expressive 3D map representation methods than 2D-OGM have been proposed for autonomous vehicle navigation [Triebel et al., 2006, Souza et al., 2013, Souza and Gonçalves, 2016, Pfaff et al., 2007, Steinbrücker et al., 2014, Hornung et al., 2013, Yue et al., 2016, Funk et al., 2021, and these methods can be broadly categorized as elevation-maps or full 3D pointcloud maps. The map's scale in robotics society is not necessarily very large considering that robots only operate in an area of relatively limited size.
In the autonomous driving society, unlike the map representations in robotics field, HD-maps can be extremely large, which can be up to a city-scale. Therefore, high efficient and compactable map representations are desirable for autonomous vehicle navigation. The recent road-based mapping methods for autonomous driving in urban environment [Schreiber et al., 2013, Ranganathan et al., 2013, Lu et al., 2017, Rehder and Albrecht, 2015, Jeong et al., 2017, Qin et al., 2020, Herb et al., 2019, Qin et al., 2021 belong to this line of research, where varieties of semantic 2D and 3D features are segmented from road region. The 2D elements are basically on road surface, including lane markers and curb-lines, zebra-crossings, and all kinds of ground arrow signs etc. The 3D structures are mostly regular vertical objects, such as traffic lights, traffic signs and light poles etc. Generally, deep learning methods are used to segment these semantic image regions.
Then, high-precision GPS/IMU and multi-ray LiDAR sensors are combined to build a prior urban map. Thereafter, the locations of the extracted semantic-feature in images are transformed into the world coordinates and anchored to the created prior map. Given the created HD map, one can apply it to accurately localize a vehicle with a single camera (aided by a low-cost GPS receiver and wheel-odometers for coarse localization in general), where the image from the camera is segmented by a deep-learning model to find the road-scene semantic features, and thereafter a fine localization is achieved by feature association and PnP (perspective-n-point) algorithms.
By comparing the above two types of HD-maps in robotics domain and autonomous-driving industries, respectively, we find that neither of them can well meet the requirement for creating HD-maps of general urban environment, where semantic road markers or regular man-made road-scene infrastructures (e.g., traffic light poles) might be absent, and there might also exist multiple-level road surface (e.g., underpass and multiple overpass). In robotics field, the elevation-based methods cannot fully represent complex scenes and the full 3D point-cloud based maps generally consume too much memory or storage space.In autonomous driving society, the existing road-based mapping methods heavily rely on semantic road features. Although these features are salient and common in high-way road scenes, they are not necessarily present universally in general urban environment, e.g., various residential communities, university campuses or industrial parks etc. Thus, one has to collect new images of semantic road-scene features which are specific to each individual scene before learning a new semantic segmentation deep learning model, for which it is a heavy workload. In addition, these methods cannot guarantee to build maps in an online mode, and some of them also need expensive high-precision GPS/IMU sensors.
Given the problems of the above map representations, it is necessary to think of an efficient map which is parsimonious in memory consumption and hard-drive storage, and highly descriptive for complex environment as well. Summing up all the above thoughts, the main contributions of our mapping work can be as follows.
(1) We propose an efficient online mapping method which can be applied to general 3D urban environments to create an HD-map for autonomous vehicle navigation. The result of mapping is a parsimonious map representation by embedding sparse local point-cloud encodings and local-OGM induced traversable-region labeling.
(2) We can handle the challenging problem of creating local 2D-OGM in non-structured urban scenes based on a real-time delimitation of traversable and curb regions in LiDAR point cloud. (3) We can achieve accurate 3D mapping in multiple-layer urban road scenarios and solve the map update problem of the existing multiple-level elevation-maps caused by dynamic objects present in the mapping process by a probabilistic fusion scheme.
We have compared the proposed map representation with the popular OctoMap map representation [Hornung et al., 2013]. Empirically, the resulting map by the proposed map representation is much smaller than that of the OctoMap representation. Based on the runtime memory coverage, our map representation can afford the online large-scale map creation while the OctoMap cannot. We have also conducted extensive experiment on online map update and map-based localization, by which we have shown the accuracy and efficiency of the proposed map representation.
In contrast to the existing mapping methods, our method belongs to the more broad category of traversableregion based mapping algorithms, where the traversable region is not limited to structured urban road where the man-made semantic road elements (road markings/signs, and traffic lights etc.) are available. Actually, our method can be applied to more general natural scenes instead of just urban environment. Our method only uses a multiple-channel LiDAR sensor (VLP-16) to generate the LiDAR Road-Atlas. Noted that we do not use cameras in this work although we occasionally show images in the paper, and the purpose is to tell what the evaluation scenarios look like. Because our mapping is an online process, it can also be easily extended to some robotic exploratory tasks in an unknown environment.  [Triebel et al., 2006]). Refer to our introduction part for more details of these methods.

Related works
Representation of traversable region -In representing traversable regions, some researches [Rehder and Albrecht, 2015, Jeong et al., 2017, Qin et al., 2020, Herb et al., 2019, Qin et al., 2021 focused on the process of building road maps. Regder et al. [Rehder and Albrecht, 2015] detected lanes in image and exploit odometry information to generate local grid maps, and pose between local grid maps is further optimized and stitched. Jeong et al. [Jeong et al., 2017] first differentiated road markings, and more informative class information can be used to avoid ambiguity on matching. Pose graph optimization based on loop closure is performed to eliminate drift for global pose consistency. Qin et al. [Qin et al., 2020] utilized road markers to build a semantic map for underground parking lots with application to autonomous parking based on only cameras. Herb et al. [Herb et al., 2019] proposed a crowd-sourced solution to build a semantic map. However, it was difficult to be applied because its computational complexity is generally high in matching inter-session features. Very recently, Qin et al. [Qin et al., 2021] proposed a light-weight semantic map representation for assisting camera-based vehicle localization in high-way structured environment. In utilizing the created road HD-map, Schreiber et al. [Schreiber et al., 2013] proposed to detect curbs and lanes, and then match the structure of these elements with an HD-map. Ranganathan et al. [Ranganathan et al., 2013] utilized detected corners on road markers for localization. Lu et al. [Lu et al., 2017] formulated a non-linear optimization to match road markers with the map. The vehicle odometry and epipolar geometry were also taken into consideration to estimate the 6-DoF camera pose. Morales et al. [Morales et al., 2009] proposed an autonomous system for robot navigation in outdoor cluttered pedestrian walkways. As a part of this system, a module of driving path detection is given for safe navigation of robot. Different from our method, this method applied a cascaded height filter to find traversable region in each individual scan, and then detected curb edge point to construct a map for navigation. The navigation map is still two-dimensional because the planned way-points lie in a planar map.
Representation of 3D space occupancy -Occupancy elevation map (OEM) is also a frequently used map representation method to deal with the problem of under-representing 3D environment in 2D-OGM [Bares et al., 1989, Triebel et al., 2006, Kweon et al., 1989, Souza et al., 2013, Souza and Gonçalves, 2016, as shown in Fig.2. In OEM, the height of the surface is stored in the corresponding cell of the discretized grid. While the OEM provides an efficient map representation, it lacks the ability to fully represent vertical structures or even multiple levels of horizontal planes (e.g., underpass/overpass bridges in Fig.2(b)). The extended occupancy elevation map (E-OEM) [Pfaff and Burgard, 2005] can alleviate the problem of OEM to some extent, and allows to handle overhanging objects (e.g., tree branches) or represent the underpass. But it cannot fully represent the environments where there are more than one level of horizontal surface (e.g., the overpass bridges in Fig.2(c) or multiple-level parking garages). To address this issue, Triebel et al. propose multi-level surface occupancy elevation map (MLS-OEM) data structures [Triebel et al., 2006], where a list of surfaces are stored in each cell of the discrete grid (e.g., Fig.2(d)). Additionally, they use intervals to represent vertical structures. The MLS-OEM roughly provides a 3D representation without the complexity of a full 3D point-cloud map, enabling a mobile robot to model complex environments. However, the MLS-OEM representation is designed for mapping static environment, and it is not capable of handling the influence of dynamic objects during the mapping process. In constrast, our method can not only represent the multiple-level surface elevation but also deal with the influence of dynamic objects in the mapping process.
Another widely used map representation is the OctoMap [Hornung et al., 2013], where each point can be stored and indexed in an Octree data structure. The indexing of each point is fast due to the tree structure. The OctoMap map is able to model almost any kinds of environments without prior assumptions on it, and it models occupied areas as well as free space. Unknown areas of the environment can also be implicitly encoded in the map. Although OctoMap has high representative capability, its memory usage and hard-drive storage is also quite high. The normal distribution transform traversability mapping (NDT-TM) representation [Ahtiainen et al., 2017] for outdoor environments is able to classify sparsely vegetated areas as traversable without compromising accuracy on other terrain types. The NDT-TM representation is proposed by exploiting 3D LIDAR sensor data to incrementally expand normal distributions transform occupancy (NDT-OM) maps [Saarinen et al., 2013], which is a 3D spatial model based on a regular grid that concurrently estimates both the occupancy and the shape distribution in each cell. In addition to geometrical information, the NDT-OM representation is augmented with statistical data of the permeability and reflectivity of each cell of NDT-OM map. Using these additional features, a support-vector machine classifier is trained to discriminate between traversable and nondrivable areas of the NDT-TM maps. Despite of advantages over 2D grid mapping, the construction of full 3D occupancy grid models [Hornung et al., 2013,Ahtiainen et al., 2017,Saarinen et al., 2013 typically has excessively high demands for run-time memory and hard-drive storage during its online application to a mobile robot in large-scale outdoor scenarios or when fine map resolution is necessary. The remainder of the paper is organized as follows. We introduce our method in Section III, and the implementation details and experimental results are given in Section IV. We discuss our method in terms of its scalable, extendable, and applicable properties in Section V. The conclusions are drawn in the last section.

Our method
Our work aims at building a 3D global sparse occupancy grid map with labeled traversable regions and static obstacles (e.g., trees and buildings) codings embedded inside. To do that, we estimate the LiDAR's 6D pose based on a LiDAR-based SLAM or odometry method. At each LiDAR's key-frame location, we build a local occupancy grid map and align it to the estimated 6D pose (Fig.8). A stitching operation is given to fuse two overlapping local occupancy grid maps when necessary. We consider not only common outdoor scenes but also the scenarios where it is necessary to build multiple layers of local occupancy grid map, e.g., the crossovers and multiple-layer parking lots shown in Fig.1.
Our method takes in point cloud data from a 16-channel LiDAR sensor and estimates relative pose between two consecutive frames. An overview of our method is shown in Figure 3, where the point cloud is processed and automatically labeled, with each point being classified as traversable regions or obstacles. Thereafter, the processed point cloud is used to construct a local grid map. Meanwhile, for each frame of point cloud, spare descriptors of local geometric structures are extracted and embeded into the global 3D LiDAR Road-Atlas along with the labeled traversable region in each local grid map.

Detecting traversable area and obstacles in each LiDAR point-cloud keyframe
Accurate detection of traversable regions is crucial for path-planning and route-following in the "repeat" stage of the "teach-and-repeat" navigation paradigm. Especially when localization in the map is not robust enough, a reliable traversable road extraction method can be able to keep the robot stay and run in the path. Traversable region labelling is strongly relevant to detection of obstacles, ranging from nontrivial ones such as vehicles, walls, overhanging tree branches and overpass bridges to trivial ones such as curbs or ditches.
Our algorithm is designed to automatically detect traversable region during map building process. The online detection is also valuable to planning paths for robotics exploration task. In general, one usually uses the Random Sampling and Consensus (RANSAC) algorithm [Fischler and Bolles, 1981] to find a dominant plane of the scene. However, there usually exist multiple planes of similar size in a scene. When applying a single RANSAC, small obstacles (e.g., curbs or pits on road) are usually overlooked if a large tolerance threshold is adopted. As a result, it is neither reasonable nor robust to apply a single RANSAC model to detect traversable regions.
Although a single RANSAC does not work well in a whole point-cloud frame, it can do a good job locally. Therefore, one could apply multiple RANSAC models sequentially to the whole point cloud, where a partial traversable region is detected and removed from the point cloud by a different RANSAC model in an iterative way. Although this multiple-RANSAC way could generally get better results than a global RANSAC model, it is not trivial to determine the exact number of iterations necessary for a good segmentation.
We propose a new multiple-RANSAC based method to detect traversable and obstacle region from a frame of point cloud. Specifically, we divide the point cloud P into overlapping local sectors A i , and perform multiple RANSAC models on these sectors, respectively, to separate traversable areas from obstacle regions. Note that the smallest obstacles by our method can range from curbs to small ditches or stones on road. Figure 4 illustrates the way of splitting a whole LiDAR point-cloud frame into many overlapping regions. Based on the splitting, the multiple-RANSAC traversable region and obstacle detection method is given in the Algorithm 1.
Since RANSAC cannot distinguish between the ground and the wall, we define check normal to check the normal vector of the plane obtained by RANSAC to remove too steep planes. The return value of check normal is only true when the angle between the normal vector and the z-axis is less than the threshold. The threshold depends on the vehicle. We set it to 0.4rad in our experiments.
Generally, there are two possibilities in fitting a plane in a split local region. One is that the whole local region can be perfectly fitted by a plane model, e.g., the left dark rectangle shown in Fig.5(a)(b). The other is that partial of the split local area can be fitted by a plane model, e.g., the red rectangle and the right dark rectangle shown in Fig.5(b). For such local regions, the plane is detected as the part which covers more than 50% of the whole local region. In Fig.5(a), we show an alternative local-region splitting scheme where there is no overlap between two neighboring local region, for which the segmentation of small obstacle is worse than that in the splitting case with overlap. Therefore, the local-region splitting way with overlap is our Algorithm 1: Detecting traversable area and obstacles Input: whole point cloud P , local point cloud choice. It deserves pointing out that the splitting with overlap is crucial for reliably detecting small objects such as curbs, which can guarantee creating an accurate local occupancy grid map, as shown in Fig.6.
In our implementation, we split a whole point cloud frame into local areas using the LiDAR image representation . The raw output data of a multiple-channel LiDAR are actually based on spherical coordinate system, including the yaw angle θ of each laser beam and measurement distance d. The pitch angle φ of each laser beam can be loaded from the calibration file. Actually, there are a couple of ways to divide a point cloud to meet the overlap condition. However, the raw point cloud representation is sparse, and it is not straightforward to directly perform neighborhood calculation and extract the neighborhood that meets the overlap condition. That is why we choose the LiDAR image representation in traversable region detection.  Figure 5: A demonstration of the role of our "'splitting with overlap" scheme on segmentation of a curb-like region. (a)Two split regions without "overlap" between them. (b) Three split regions with "overlap" between two neighboring local regions. Based on Algorithm 1, the total traversable region in a point-cloud frame is composed of many smaller traversable planes which are detected in the local regions. The rest corresponds to obstacle area. In (a), part of road plane (blue region) is detected as the curb obstacle. In (b), we can deal with this issue.

Creating 2D local occupancy grid map from each LiDAR point-cloud keyframe
Conventionally, the 2D-OGM is obtained by a 2D LiDAR sensor. It can produce a clear boundary representing dynamic or static obstacles in indoor environment or 2D structured outdoor environment, and thus a 2D-OGM can well represent the environment. In such environment, vehicle localization and path planning can be fulfilled with ease in a 2D-OGM generally. In a more general 3D urban environment, 2D-OGM cannot be reliably created if obstacles are not well labeled. In our work, we can deal with this issue by the traversable and curb region labeling procedure as described in the last section. Meantime, we can also take the advantage of 2D-OGM over its 3D counterpart, i.e., efficiency and small memory usage.
Specifically, we assume that the local ground area around the vehicle can be approximated as a 2D plane. After removing points with significant height, each frame of point cloud is projected vertically onto this 2D plane. Figure 6 (a) shows such a projection. However, the resulting points are too dirty to apply the 2D occupancy grid mapping algorithm [Thrun et al., 2005, Hess et al., 2016 to construct a local occupancy grid effectively. Therefore, we hope to convert all the projected 2D points to obtain a virtual scan which is more convenient for generating a 2D-OGM before applying the method in [Thrun et al., 2005, Hess et al., 2016.
To do that, we adopt the ray-tracing technique to find the closest obstacle point along each ray direction, e.g., the yellow points in Fig.6(a). In addition, if there is no obstacle point along a certain ray, the farthest point, which has been labeled as traverable (e.g., the blue points in Fig.6(b)), is used to speed up the filling of the map to differentiate known area from unknown one. Figure 6 shows a created 2D-OGM. In our implementation, if the resolution of our 2D-OGM is set to 0.1m/cell-width and the size of the local region that LiDAR can sense is 40mx40m, the size of our 2D-OGM is 400 pixels by 400 pixels.
Due to measurement noise, it can be assumed that the altitude of each local grid of the segmented traversable region in each LiDAR frame is subject to a Gaussian distribution. Specifically, the mean of each Gaussian distribution is the height of the local grid that is represented in the global coordinate system. In general, the standard deviation of each Gaussian distribution increases with the distance, d, to the LiDAR sensor. For efficiency purpose, we simply set the standard deviation as a combination of two parts, namely σ = α d + β.
The α d is a variable standard deviation which becomes larger with the increased distance d, and β is a constant used to ensure that the standard deviation of the cell under the robot is nonzero. Figure 7 shows an example of the distribution.

Building 3D global occupancy grid map
With a LiDAR-based odometry method, we can obtain the pose trajectory of the vehicle. For each key-frame pose along the trajectory, we estimate a local 2D-OGM and make it aligned with the pose represented in the global coordinate system. To construct a 3D global occupancy grid, our solution is to merge all local occupancy grid maps by stitching the neighboring ones so that the transition between two local 2D-OGM is as smooth as possible, as shown in Fig.8.
To make our explanation easier to understand, let us use Fig.9 as a conceptual example although the real scenarios should be different from it. As shown in Fig.9, the scene includes an underpass and an overpass. When a robot carrying a LiDAR sensor moves along the dark line in Fig.9(a), correspondingly, the whole route can be split into five segments (red→dark→green→dark→blue, shown in (b)) based on the visibility of the cell m to the LiDAR sensor. For example, the red, green and blue segments are the first, third and fifth one, respectively, where the cell m is visible to the LiDAR sensor. The two dark segments represent the two courses where the cell m is non-visible to the LiDAR sensor, and should not be considered for fusion. Note that there are two levels at the cell m location, and the top level of the cell m is visible only to the red and blue segments, and the bottom level of m is visible only to the green segment. To efficiently generate the 3D-OGM, only the representative observations of the red, green and blue segments are utilized during fusion. Let us take the representative observations of the red segment as an example. When the robot moves  from the starting point to the end of the red segment, the distance between the robot and the cell m first decreases and then increases, with the minimal distance achieved at location A. Generally, the variance of measuring the ground cell m at location A is the smallest. We set the observation made to the cell m at location A is denoted as the representative one of the red segment. Likewise, the observations made to the cell m at locations B and C are the representative ones of the green and blue segments, respectively.
We first give some notions before the detailed description. A certain cell in the real world (the global map) can be denoted by c i , and there might be multiple-level planes at c i . The distance between the moving robot and c i over time is calculated as d ci (t), where t is a time variable. The observations made to c i by a robot over time in a segment s j (e.g., the red one in Fig.9) are denoted as m sj ci (t), t ∈ [the period when robot was running in s j ], where the representative observation is represented by m sj ci (t * ). That means the robot is closest to c i in the segment s j at time t * .
As illustrated in Fig.9, if assuming the whole course includes N segments, s j , j ∈ [1, N ], where the cell c i is totally visible to the LiDAR sensor, we thus have N representative observations made to c i , m sj ci (t * ), j ∈ Before we estimate the GMM representation of the altitude at each cell c i , we would first apply a merging step to remove the redundant ground planes. For example, the representative observations made at locations A and C to the cell m in Fig.9 should be merged because they are at the same altitude. In contrast, the representative observations made at A and B, or C and B, respectively, are not at the same altitude and should not be merged.
According to the method in [Sun and Wang, 2011], we use the overlaprate (OLR) to evaluate whether these independent observations should be merged or not. To do this, the set of 2D local occupancy grids that correspond to the representative observations m sj ci (t * ), j ∈ [1, N ] are used, represented by ogm sj ci (t * ), j ∈ [1, N ]. Note that the altitude of each of these occupancy grids can also be represented by a Gaussian mixture model in the global coordinate system (refer to Section III.B). Thus, we first sort them based on the means of their altitudes.
The lowest one with the smallest mean value of altitude, e.g., ogm sq ci (t * ), is initially assigned the label "1". The second lowest, e.g., ogm s k ci (t * ), is compared with ogm sq ci (t * ) based on altitude difference, and is assigned a label by checking whether the altitude difference is larger than a preset threshold. Generally, given the label of ogm sq ci (t * ), f sq ci , the label of its adjacent neighbor, f s k ci , can be updated based on (1), where OLR(·, ·) = p saddle /p sub max . The p saddle is the saddle-point location of the joint distribution of the altitudes of ogm sq ci (t * ) and ogm s k ci (t * ), and p sub max is the location of the smaller peak of the two modes. The is an threshold and can be 0.6 as described in [Sun and Wang, 2011].
In the end, we get Z independent labels which means there are Z layers of ground surface at c i . Then each observation through the whole course can be assigned one of the Z labels based on (2).
where µ t is the altitude of the observation made at time t, and the variables with * are associated with the representative observations. For example, the t * is the time instance at which a representative observation is made, and µ * is the mean altitude value of a certain representative occupancy grid.
Finally, we will get Z sequences that can be used for fusion. According to the general elevation map fusion method [Triebel et al., 2006], the Z sequences are fused sequentially to generate the final result. Note that because our framework is an online mapping process, it is not necessary to wait until all observations thought the whole course are made. The update in (3) is an online process as well in real implementation.
Where µ f,0:n f and σ f,0:n f represent the fusion result from the first observation to the last observation that has been assigned the label f , and n f is the number of observations labeled as f .

Embedding geometric information for vehicle localization
The above steps can be used to create a 3D road-map or traversable-region map, which is mainly for vehicle path planning purposes. Sometimes it is also useful to keep vehicle run in the track when vehicle localization is not accurate or totally gets lost. But for accurate vehicle localization, the above created map cannot be reliable generally although it should be applicable to 2D structured scenes. When vehicle is running in a non-structured 3D scene, we have to incorporate 3D vertical geometric structures into the created road map. Inserting vertical structures into the road map is in analogy to growing trees along road which stretches away into the distance.
To do that, vertical structures should be detected in each LiDAR point-cloud frame and then inserted into the generated road-map. Based on our scheme presented in Section 3.1, we initially get the pointcloud region pc obst which corresponds to obstacles (e.g., trees, bushes, light poles or buildings). However, Dynamic objects (moving pedestrian or vehicles) can also be detected as obstacles because this scheme cannot distinguish dynamic objects from static obstacles. Dynamic objects are not expected to be reserved as vertical structures in the final LiDAR-Atlas because they can affect adversely localization performance when the LiDAR-Atlas is used for navigation. In addition, the static obstacle points are somewhat redundant for localization and should be sparsified for efficiency purposes.
To deal with these issues on vertical geometric structures, we propose a probabilistic way to achieve the aim of reducing redundant 3D points and removing dynamic objects as much as possible. Specifically, an efficient 3D probabilistic grid map, 3D-PGM, in x−y −z dimensions (z representing the vertical dimension) is created to reserve sparse vertical static structures, where the resolution of x − y plane is set to be 0.1m/cell-width (the same as the one adopted in creating the 2D-OGM in Section 3.2). Vehicles actually drive on the ground, and there is no significant change in height for dynamic objects in vertical dimension locally. So we set a sparser resolution in the vertical dimension. The lowest grid in the 3D-PGM is slightly higher than the ground so as to effectively ignore the noisy obstacle points that may exist near the ground.
To remove dynamic obstacles through the 3D probabilistic grid map, we find that it is not necessary to use high-precision numbers to represent probability values. Specifically, for each 3D probability cell, we can use fewer bits (4 bits) to store a probability value to reduce memory cost during runtime, shown in Fig.11 (a). The encoding process is a simple linear mapping of probability values to the range 1 to 15. Exclusion of 0 is to ensure that there are an odd number of states, ensuring that 0.5 (unknown state) can be represented exactly. So 1 to 7 means non-obstacle, 8 meaning unknown, and 9 to 15 meaning obstacle. In contrast, 2-bit representation cannot express uncertainty under the condition of odd state (e.g., the part greater than 0.5 or less than 0.5 contains multiple possible values), and 3-bit representation cannot fill all bytes when the number of segments is not a multiple of 8. 4-bit representation is the most saving and appropriate option. Then, we encode the occupancy probability for all segments at the same 2D position based on the 4-bit representation. In implementation, the probability update is by mapping back to the range of 0 to 1 from the 4-bit representation based on efficient look-up table. Considering the result of this encoding as a descriptor, the representation is ultimately consistent with a 2D map.
(a) (b) Figure 11: Grey-scale level shows the occupancy probability of a segment. The darker the grey-scale level, the larger the probability of being occupied. (a) is an illustration of probabilistic coding, where multiple segments corresponding to the same 2D location are shown. In this example, we encode the probability with 4 bits and label it on the right. (b) shows the process of converting an obstacle point cloud to a 3D probability grid map for a LiDAR frame. The dark dot is the position of the LiDAR sensor and the three red dots represent three obstacle-grid locations. Once we have the economical representation of probability values, given a LiDAR frame, the point cloud pc obst (e.g., the red region in the Fig.10 (b)) corresponding to the vertical structures (including dynamic objects) can be represented with the 3D probability grid map. This procedure is similar to the one used (a) Before the car leaves.
(b) After the car leaves. Figure 13: Illustrating dynamic object filtering by probabilistic fusion: A car is leaving the parking-lot.
for constructing a general 2D occupancy grid map, where a miss probability p miss is applied to all segments passed by the ray from the LiDAR to the obstacle point, and a hit probability p hit is applied to the segment where a obstacle point is located. For ease of calculation, the sensor model we use takes fixed values for both p miss and p hit . The process of adding obstacle points within a single LiDAR frame is shown in Fig.11 (b).
Reversely, we can decode the probability of each segment from the descriptors of each 2D cell, and based on these probabilities, we can approximately reconstruct the obstacle point cloud with the obstacle segments.
The reconstruction result of one frame is shown in Fig.10 (c). Over time, we can construct a global map from each consecutive local frames by applying a probabilistic mixture model to fuse two 3D probabilistic grid maps. Fig.12 shows a update of the probabilistic grid map by adding a new frame (each key frame as a local map) to the existing global map. We also show one example of filtering out dynamic objects with the probabilistic mixture in Fig.13.
To apply the map to robot localization, one can reconstruct a point cloud from the 4-bit descriptors of a LiDAR frame, and utilize ICP-like registration method for vehicle localization in the map. In our localization implementation, we first encode the frame to be localized, then decode the point cloud for ICP, and then do ICP matching with the point cloud decoded from the global map. The purpose of encoding and then decoding the LiDAR frame to be localized is to make it consistent with the sparsity of the global point cloud. Figure 14: Examples of creating local 2D occupancy grid maps (2D-OGM) in real environment test, with the columns corresponding to the images captured by a camera installed on the robot, the frames of point cloud by VLP-16, the detected curb regions, and the created 2D-OGMs based on the curbs, respectively.
We adopt the VLP-16 (16-channel) LiDAR as the only sensor to implement the proposed method, although we believe more advanced multi-channel LiDAR sensors (e.g., 32 or 64 channels) should work equivalently or even can achieve better results.

Traversable-and curb-region detection accuracy
In general, our method is a road-based or traversable-region based map representation. Therefore, traversable-and curb-region detection performance is first evaluated on the KITTI-360 benchmark dataset [Xie et al., 2016]. The reason we choose the KITTI-360 dataset is that this benchmark dataset provides ground-truth semantic segmentation of street scenes. Curb detection is crucial for autonomous driving in urban scenes. Actually, in most urban areas, the lowest positive obstacles no matter for autonomous cars or for robots are probably curbs, where positive obstacles are defined as the objects which stick out above the ground. The other positive obstacles such as vehicles or pedestrians are relatively easier to be detected. In addition, we also need curb-region detection results to build the 2D local grid map (Fig.6).
Therefore, both traversable-and curb-region detection performance is evaluated to ensure the accuracy of the created road map. Although our framework is implemented with a VLP-16 LiDAR sensor, it should also work with the other multi-channel LiDAR sensors. To evaluate with LiDAR sensors of different number of channels, we tried our method on the original point cloud data of the KITTI-360 benchmark and on the subsampled versions, respectively. Specifically, we sample one channel from every four channels of the KITTI-360 LIDAR data to simulate a 16-channel LiDAR sensor. We sample one channel from every two channels of the KITTI-360 LIDAR data to simulate a 32-channel LiDAR sensor. The KITTI-360 benchmark provides multi-frame 3D semantic labels. We divide the point cloud into cells at the resolution of 0.1mx0.1m per cell, and label these cells as road-, curb-, or irrelevant-cells based on semantic information, where the irrelevant cells mostly contains non-trivial obstacles such as vehicles or trees. Then our curb and traversable detection results represented in the 2D cells are compared with the ground-truth.
For curb region, we evaluate the N missed only, which is calculated as the number of curb cells that has been labeled as traversable. Additionally, we deem it as correct if a curb cell is labeled as an irrelevant-cell. For traversable region, we evaluate the N f alse only, which is the number of road cells that are mislabeled as obstacle. Given the total number of road cells N road and total number of curb cells N curb , we define the missed detection rate of curb cells rate missed = N missed N curb and false detection rate of road cells rate f alse = N f alse N road . Based on different threshold in RANSAC, we finally get the performance curves as shown in Figure 15. From the figure, we observe that the missed-detection rate deceases with more channels of LiDAR beams. Counter-intuitively, we have noticed that the false detection rate of traversable-region cells increases when using the LiDAR sensor of more channels. The main reason is that the LiDAR sensor with more channels has more points, so there are more false road points generated in general, and these points will fall into more cells. As a result, more cells are incorrectly marked as curb.

The 3D LiDAR Road-Atlas Results
To evaluate the whole 3D mapping process, we have evaluated our approach in three real-world urban scenes and two virtual urban scenes (garage and campus). The three real-world urban scenes include our own campus scene and two scenes from the KITTI Odometry datasets. The two virtual urban scenes were parts of the Autonomous Exploration Development Environment [Cao et al., ] built by a CMU robotics team. All evaluations are based on a 16-channel LiDAR sensor. Figure 16: The obtained global 3D LiDAR-Atlas of five scenes. The top three rows: the first column shows a picture of the whole scene except the two KITTI odometry scenes, for which we do not have the global pictures. The second column shows two views of the created LiDAR road atlas for each scene. To perceive the altitude variation at each location of the created 3D road map, we have encoded the altitude of each cell in gradually changing colors, where red color represents the lowest in altitude, blue representing the highest in altitude, and green and yellow in the middle. The third column shows this effect. The last row shows correspondingly the LiDAR-Atlas resuts for the two KIITI-odometry scenes.
In the three real-world urban scenes, we use the LeGO-LOAM approach [Shan and Englot, 2018] to estimate a relative pose between two consecutive LiDAR frames, and use the LiDAR-Iris method [Wang et al., 2020] to provide a stable loop-closure detection. In the two virtual scenes, we use the ground-truth pose estimation Figure 17: The embedded static-obstacle points for vehicle localization (red) and the traversable region for vehicle navigation (white). The scenes are: our campus, v-campus, v-garage, KITTI-00, KITTI-05.
provided by the datasets. In Fig.16 and Fig.17, we show the obtained 3D LiDAR road atlas without embedded vertical structures and the 3D LiDAR road atlas with embedded vertical structures, respectively.
Note that although the tested real-world scenes contain non-flat road scenarios, there are no multiple-layer road infrastructures, such as underpass and overpass. In contrast, the simulated garage and campus scenes provided by the Autonomous Exploration Development Environment [Cao et al., ] contain such multiplelayer road or traversable-region scenes. The obtained LiDAR Road-Atlas maps from the two simulation environments have shown that our method is effective in dealing with such scenes. In implementation of traversable and curb detection in the simulation environment, due to the large horizontal resolution of the LiDAR in the simulation environment (350 points per circle), we adjusted the width of LiDAR imagery to 350, and the width of the local window for detecting traversable areas and obstacles to 9 (also about 10 • ).

Memory usage and map size
The total memory consumption of our method includes two major parts, i.e., the cost for the created 3D road map (e.g., the results shown in Fig.16) and the space required for saving the encoded vertical structure of each cell. As in the occupancy elevation map [Souza and Gonçalves, 2016], our 3D road map also takes the advantage of the probabilistic fusion when merging the 2D local occupancy grid maps. Generally, it can save a lot of memory compared to the other 3D occupancy grid maps, e.g., the OctoMap method [Hornung et al., 2013]. To make a comparison with the OctoMap method in terms of memory and hard-disk usage, we have recorded the file size and memory usage by our method and the OctoMap method for different scenes, shown in Table 1. The input data of the two methods are exactly the same, and the cell or grid resolution equals to 0.1m.
In the table, the "runtime" represents the total runtime memory cost till the end of the whole online mapping process. The "H-D" represents the size of the file after the map is saved as a file in hard drive. In our implementation, the runtime memory is dependent on several major factors, i.e., the scanning range of LiDAR sensor, the cell resolution (representing the actual size in 3D world per cell), the number of vertical segments which are used for encoding the vertical structure within each cell (determining the height limit of the point cloud within each cell), the number of bits per segment and the number of collected keyframes during traveling unit distance (e.g., 1 meter in our work).
To make a fair comparison with the OctoMap representation, the cell resolution and the number of collected keyframes across each whole dataset are set to the same for both our method and the OctoMap algorithm [ Hornung et al., 2013]. Therefore, the rest parameters that can be adjusted for determining the "runtime" and "H-D" values of our method is the scanning range of LiDAR sensor and the numbers of bits used to encode the vertical structure within each cell.
In terms of the scan range of LiDAR sensor, it can affect the "runtime" memory of our method significantly. In our 3D road map, we temporarily stores the local 2D occupancy grid maps as the intermediate results in order to allow for updating poses of frames (such as the pose optimization when loop closure happens). In general, the size of a local 2D-OGM depends mostly on the LiDAR's sensing range and more memory consumption is necessary as the sensing range becomes larger. Normally this value can be set as the range limit of a LiDAR senor. For example, the VLP-16's scan-range limit can be up to 100 meters. In this case, the size of a local 2D-OGM is 2000x2000 given the cell resolution of 0.1m/cell, which means that there are about 4x10 6 cells in each local 2D-OGM.
Empirically, the point cloud 60 meters away is generally quite sparse, meaning that it is not necessary to use the whole frame of point cloud for building map. Considering these factors, it is preferable to reduce the size of local 2D-OGM by only utilizing point cloud within a certain small range. For example, we set the radius of LiDAR's scan range to be 20m and 30m, respectively, in our experiments. Additionally, we can further reduce the memory cost on the descriptor of encoding vertical structure of each cell from four bytes to one byte only. Accordingly, the runtime memory usage of our method with different settings of mapping parameters are shown in Table 1.
When only using point cloud within a scan range of 20m for building 2D-OGM, the runtime memory usage in all environments is smaller than that of OctoMap except for the virtual garage environment. This result applies no matter whether 4-byte or 1-byte descriptors are used for representing the vertical structure of each cell. We can also observe a significant reduction of runtime memory by using 1-byte descriptors in contrast to using 4-byte descriptors, where we observe a reduction of about 23%, 25%, and 17% on runtime memory cost in our campus, the v-campus and v-garage scenes, respectively. It deserves pointing out that the garage environment is a multi-layer parking-lot environment (as shown in the Figure 16) and the trajectory is more concentrated within a small space. For such kind of scenarios, the memory usage of OctoMap is smaller.
In contrast, our method needs to build the 3D road atlas by tracing a quite long trajectory within the same area repeatedly. Thus, the intermediate results (the local 2D-OGM) can occupy lots of memory, which results in a larger runtime memory consumption than the OctoMap. When using point cloud within the scan range of 30m for building 2D-OGM, the runtime memory usage of our method is generally more than that of the OctoMap method. This is still caused by the memory coverage of the intermediate results which are prepared for pose correction during loop-closure process. But we also notice that the OctoMap method will crash when applied to very large scenes, such as KITTI00 and KITTI05 odometry benchmark datasets, due to immense memory cost and extremely long map update time. In contrast, our method can still be applicable to such kind of large scenes.
In terms of the hard-drive storage space, it can be seen from Table 1 that our method saves a lot of more storage space than the OctoMap, no matter whether in the simulated environments or in real scenes. We also observe that the hard-drive storage of the map by our method does not change much when using either 4-byte or 1-byte descriptors for representing vertical structure geometry within each map cell. For all the results shown in Table 1, a 16-channel LiDAR (or simulated 16-channel LiDAR) sensor is used, where the point cloud within a radius of 30m is used for building local 2D-OGM with only 1-byte descriptor for encoding vertical structure of each cell. Because it is the completed map that will be used in the "repeat" stage of the "teach and repeat" paradigm, considering the small hard-drive storage space of our method, our map representation is more intriguing than the existing 3D LiDAR maps, such as the OctoMap, when applied to a large-scale environment.
We also compare our method with the LeGO-LOAM approach [Shan and Englot, 2018] in terms of the harddrive storage size of the resulting map file. In the LeGO-LOAM method, only points with high curvature or points within large plannar regions (corresponding to walls or road) are reserved in the resulting map. On three compared datasets, the "H-D" values of the LeGO-LOAM method is smaller than ours, but the pointcloud map obtained by the LeGO-LOAM method alone cannot be used for navigation because traversable region is not provided in the map. In contrast, the LiDAR Road-Atlas obtained by our method is a complete map which can favorably for navigation purposes alone.

Online map update
Our map representation can be applied to an online SLAM process, for which it is important to update the map instantly. Thus, we tested the time consumption on map update. In the implementation, the map is updated once per second. When there is no pose update, the update will be skipped. The frequency of the LeGO-LOAM back-end optimization for pose update is less than 1Hz. The code of the map update part is paralleled in eight threads. We record the time of map update for three scenes and get the results as shown in Figure 18.
Specifically, for the three scenes that used for evaluating the map-update time cost, the first scene is part of our campus (the first row of Fig.16), to which our full online mapping procedure is applied, including front-end pose estimation, loop-closure detection, back-end optimization, traversable-region (curb) detection, 2D-OGM generation, and multiple-layer road map fusion. The second and third scenes are the virtual campus and garage scenarios, where poses between consecutive LiDAR frames are known beforehand. Thus, we do not need the front-end pose estimation, loop-closure detection and pose optimization steps.
We can observe that most map update can be finished within one second when point-cloud data within the radius of 20m are used for creating the local 2D-OGM, shown as the bottom row of Fig.18. In contrast, it generally takes more time to update map when using point-cloud data within the radius of 30m, shown as the top row of Fig.18, which makes many senses in that a larger local 2D-OGM needs more time to be incorporated into the global map.
The first column of Fig. 18 shows the results of running SLAM online and updating the map in our campus environment. It can be seen that there are two obvious peaks near the 141 th and 181 th updates. This is due to the fact that more local map poses are updated after finding two loop closures. The other results shown in Fig. 18 corresponds to the map update in the two virtual scenes, where relative poses between two LiDAR frames are given in the simulation environment and thus no significant time-cost variations are present. The last column shows the results of the simulated garage scene. It can be clearly seen that the update time increases mostly linearly as the number of map updates increases. This is because the garage itself covers a relative small area, but it has more floor layers. The number of local maps covered in the same 2D-OGM cell continues to increase, and the time used to update a cell is proportional to the number of local maps that need to be fused.

Vehicle localization in the map
Once our map is ready, we first test the performance of vehicle localization in it given the LiDAR point cloud data. Our test of localization is performed based on iterative closest point method (ICP) by directly matching current LiDAR point cloud with the obstacle point cloud which is decoded from our map representation. Since ICP is sensitive to the initial value, we use the localization result of the previous frame as the initial value of the new frame. At the same time, in order to accurately describe the localization error, we use the ground-truth pose in mapping to exclude the pose error of the map itself caused by SLAM in the results. The way of our localization experiment is similar to that in Monte Carlo localization methods. The Range-MCL ] is a method with better experimental results in the class of Monte Carlo localization methods. It provides test results on multiple public datasets including Apollo [Lu et al., 2019] and MulRan [Kim et al., 2020] in the paper, so we compare our results with the Range-MCL. We test on Apollo and MulRan datasets, which were collected when vehicles ran more than once in the same real scenes. Therefore, both datasets are suitable for testing localization performance given the map representation. We take the same trajectories as used in Range-MCL for ease of comparison. For the Apollo dataset, since the beginning of the test trajectory is not included in the map, it is difficult to locate it effectively. We start the test trajectory on this dataset from the 200th frame. For mapping data, since dense continuous frames are not necessary for mapping, we take one frame every 20 frames. Testing is done on every frame. For both mapping and test data, we do not directly use the original 64-channel data, but simulate 16-channel LiDAR data by taking one channel from every four channels to validate the performance of localization using a 16-channel LiDAR and our LiDAR road-atlas representation. The maps constructed from the mapping data of both datasets are shown in the Fig.19, respectively.
For the test results, we obtain the translation and rotation errors, and the average time cost for each frame. In Table 2, we compare the results with those given in Range-MCL. In this test, our method adopts the following parameters: resolution=0.2m, radius=40m (because some scenarios are relatively open and small radius can result in inaccurate positioning). The obstacle point (static vertical structures) takes the range of −1m to 7m with respect to the sensor in vertical dimension, which is enough to cover most effective vertical structures. The number of segments in vertical dimension is 8, and each segment uses 4 bits to represent the probability. This means that the length of the descriptor for encoding vertical structures in each cell is 32-bit. To demonstrate the effect of different parameters on localization accuracy, Table 3 shows the effect of various number of segments in splitting the vertical dimension and Table 4 shows the effect of different LiDAR range radius and map resolutions. We observe that more segments can slightly improve the accuracy of localization, but at the same time, the time cost will increase due to the increase in the number of points.
In order to ensure that the occupied space does not increase significantly, our experiments simultaneously change the resolution and range while keeping the number of cells in the local map unchanged. Increasing the cell side length will reduce the number of points and increasing the range will bring a lesser increase in the number of points. It can be seen that the implementation shows that as these two increase, the time consumed decreases.  : Three examples of path planning based on A* algorithm in the created 3D maps by our method.

Path planning in the map
Since our map belongs to traversable-or road-based 3D map representation, where traversable-or roadregion has been explicitly represented in the map. Therefore, another function that our map can provide is the path-planning capability. Given a start location and a destination in the map, a global path-planning algorithm can be applied to find a route for vehicle navigation given accurate localization. Figure 20 shows three planned paths in three created road maps by our methods, respectively, where A* algorithm is applied to find the optimal path in each road map. It can be observed that the planned routes actually extend at different layer of the road map of the virtual campus and virtual garage scenes.

Discussion
The proposed map representation is built in an online mode by a LiDAR-based full SLAM framework. Combined with an exploratory planning algorithm, the proposed method can be used for autonomous robotic exploration in unknown environment. Certainly, the proposed method can also be used in an offline mode where the collected data are processed in a workstation to obtain the environment maps. In all, the created map representation can be scalable, extendable, and applicable in various ways.
Scalable. Our map representation can be readily for scaling up whenever necessary to expand the current map size. Since the hard-drive storage size of our map representation is small and compact enough, one preferable way to expand the map scale is to utilize the existing map that has been saved to hard-drive, and one only need to load the map file into memory, and to re-localize the current position of vehicle in the existing map, and do the same online mapping as this work. This way is better than continuously online building a large map in terms of the runtime memory usage.
Extendable Although we only use a single 16-channel LiDAR sensor to build the map, one can use LiDAR sensor with more channels, such as 32-or 64-channel LiDAR sensors. One advantage of our map representation when only using multiple-channel LiDAR is that the runtime memory consumption and hard-drive storage should approximately be invariant to the number of LiDAR's channels. Additionally, one can add heterogeneous information to the map by incorporating the other types of sensor data, such as visual features from cameras of different spectrum. To do that, one can estimate the external parameters between LiDAR sensor and camera before extracting visual feature along with corresponding point-cloud descriptor. Generally, richer information can be obtained from camera and stored in the map. Thus it can result in more accurate vehicle localization if vehicles are equipped with cameras as well.
Applicable Our map representation can also be used for detecting dynamic objects given the map, and also be useful for local motion planning. Because our map's hard-drive storage is small, it can be applied to distributed multiple-robot online mapping, and robots can share already saved maps with each other. Additionally, since our map representation is compactable in terms of the cell resolution, where robots can use a coarse-to-fine manner to exchange intermediate map results under bandwidth limitation whenever necessary.

Conclusion
We propose a compactable and efficient 3D map representation for autonomous robot or vehicle navigation, called the LiDAR Road-Atlas. It can be generated by an online mapping framework based on incrementally merging 2D local occupancy grid maps, where automatically labeled traversable regions, and 3D landmarks are encoded and embedded into the LiDAR Road-Atlas. Given the LiDAR Road-Atlas, one can achieve accurate vehicle localization in a coarse-to-fine manner, path planning in the labeled traversable regions and some other tasks. We compare our map representation with a couple of popular map representation methods in robotics and autonomous driving societies, and our map representation is more favorable in terms of efficiency, scalablity and compactness.