3D Environment Mapping with a Variable Resolution NDT Method

Abstract

With the continuous development of the 3D LiDAR (Light Detection And Ranging) mapping algorithm and its application in various fields, the size of the point cloud map becomes a bottleneck that limits the 3D LiDAR mapping algorithm from running for a long time. In this paper, a 3D LiDAR mapping method based on scan-to-map and variable resolution NDT (normal-distributions transform) registration is proposed. When updating the global map, variable resolution processing can reduce the size of the global map and improve the accuracy of map construction. In addition, the size of the map created by the mapping algorithm is proportional to the size of the space and does not grow infinitely over time. The mapping experiments using a rotating LiDAR in the room, corridor, and outdoor environments show that the algorithm has higher mapping accuracy and smaller map size than without considering the variable resolution strategy. The experimental results of the map construction for a long time in an appropriate test area illustrate that the map built by the algorithm does not grow infinitely with time in the fixed space. In summary, by adjusting the map resolution adaptively according to the curvature of different areas in the 3D LiDAR mapping process, the proposed variable resolution strategy can maintain the size of the global map almost proportional to the size of the space. Moreover, the mapping accuracy can be improved as well.

1. Introduction

3D environment map construction technology has a wide range of application prospects, such as precision agriculture [1], building reconstruction [2], search and rescue [3], tree surveying [4], terrain mapping [5], historical building appearance mapping [6], and so on. Therefore, it is of great significance to study the 3D environment map construction technology.

LiDAR (Light Detection And Ranging) and camera are mainly used for 3D mapping. For example, a monocular camera is used for real-time UAV (Unmanned Aerial Vehicle) dense mapping [7], an RGB-D camera is used for indoor unbounded RGB-D dense point cloud 3D reconstruction with high accuracy [8], and a LiDAR is used for environment map construction to cope with challenging agricultural scenarios [9]. The camera can only be used in a well-lit environment, while the use of LiDAR is almost unaffected by the environment. Compared with the camera, LiDAR has higher accuracy and a larger field of view, so we use 3D LiDAR to conduct research on 3D map construction.

3D environment maps are represented in many different ways, including point cloud [10], mesh [11], surfel [4], grid [12,13], polarized LiDAR Map [14], etc. The point cloud is an effective way to represent the real environment. With the rapid development of LiDAR, there is more and more research on LiDAR-based mapping algorithms, such as [15,16,17]; all are based on research of 3D LiDAR mapping algorithms. Since the point cloud can represent the surrounding environment simply and intuitively, it is more convenient to use than other types of map representation, so we use point cloud as the map representation.

When using the point cloud as the map representation, the map will become denser and denser over time, even with the limited environment. The conventional methods to solve this problem include removing redundant points according to the distance between point clouds [18,19] or the time interval over the capturing time of the point cloud [20]. Some researchers also tackle the above question by detecting the overlap degree of sub-maps and carrying out map fusion [21,22,23]. The method of removing redundant points in the map according to the distance between point clouds is similar to uniform downsampling, which is convenient to implement, but it cannot distinguish objects with different features. The strategy of removing redundant points according to the time when the point cloud is added to the global map does not make use of the previous point cloud information, which is not suitable for constructing a point cloud map. The method of detecting the overlap degree of the sub-maps requires maintaining the sub-maps, which is difficult to implement.

To solve the problem of significant redundancy caused by repeated scanning of the same space, we propose to process the map with a variable resolution strategy and adaptively adjust the point cloud density of different locations in the map according to their curvature. Thus, different locations in the map have reasonable point cloud density. This ensures that the map size varies with the size of the space rather than growing infinitely over time. In addition, in the process of variable resolution processing of the map, redundant points are removed according to the probability density value of each point in the global map, which is helpful in improving the accuracy of map construction.

The main contributions of this paper are as follows:

• A novel strategy of variable resolution processing for point cloud maps is proposed, which can reduce the map size effectively without decreasing the mapping accuracy;
• The size of the map generated with the proposed strategy increases proportionately to the environmental extent rather than the scanning time;
• The experimental results in both structured and unstructured environments demonstrated the efficiency of the proposed method.

The rest of the paper is organized as follows. Section 2 reviews some research related to this paper. Section 3 introduces the algorithm of the whole system. Section 4 introduces the experiment results and discussion. The experimental environment includes a room environment, a corridor environment, an outdoor environment, and a conference room environment. Section 5 summarizes the paper and puts forward suggestions for future modification.

2. State of the Art

This section introduces some 3D mapping systems and measures to reduce the redundant points in the map. The strategies for removing redundant points in the map include removing redundant points according to the distance between point clouds or the time when the point cloud was added to the global map, detecting the overlap degree of sub-maps and merging the sub-maps with higher overlap degree, point cloud compression, and multi-resolution maps, etc.

Removing redundant points according to the distance between point clouds is a common method to remove redundant points. There are different ways to implement it. In paper [18], before adding a new point to the global map, it is necessary to search within a certain radius of the point, and this point will be added to the map only if no results are returned. If there are already points within the radius, only the point with the least uncertainty is kept. In paper [19], the points whose distance is less than a certain threshold are directly removed. The nature of this way of removing redundant points is similar to the uniform downsampling of the point cloud, but the disadvantage is that different features cannot be distinguished.

The strategy of removing redundant points based on when the point cloud was added to the map is crude. In paper [20], a planar feature is stored in each voxel of the global map. When the uncertainty of the planar feature in a voxel converges, all the historical points are discarded, and the estimated planar parameters are retained. Once a new point appears, only the latest 10 points are kept, and the new plane features composed of these 10 points are calculated. This way of updating the map only uses the newly obtained information of a few point clouds and does not consider the information of the previous point clouds at the location. In this case, the global map is represented by plane elements, and the strategy of removing redundant points according to the time when the point cloud was added to the global map is feasible, but it may not be appropriate for the point cloud map.

It is a good way to reduce the size of the map by merging sub-maps. The map update strategy in paper [21] is to search all the sub-maps involved in the closed loop when a closed loop is detected, and then fuse these sub-maps together. Sub-map fusion can save memory by preventing new sub-maps from being created when the same space is revisited. This has the advantage of making the reconstruction complexity proportional to the amount of space explored, rather than the duration of the exploration. However, this scheme only works when the closed loop is detected, and searching all the sub-maps after the closed loop is detected also needs more time. Paper [22] introduces a sub-map fusion strategy based on overlapping scan space on the basis of paper [21]. Sub-map overlap is calculated by comparing volume reconstruction and occupancy information stored in each sub-map. If the overlap of two sub-maps exceeds the threshold, it indicates that there is significant redundancy between the sub-maps, and map fusion is required. The calculation of sub-map overlap also requires traversing a large number of sub-maps, resulting in a large amount of computation. In paper [23], a method to calculate the map increment is proposed, and only the incremental part is added when the map is updated. It is also a time-consuming operation to calculate the map increment in this method.

There have also been some studies on reducing the size of the map by compressing the point cloud. In paper [24], a Recurrent Neural Network is used to compress the point cloud. In paper [25], a clustering method is used to compress the point cloud. In paper [26], a voxelized point cloud compression scheme is proposed. In paper [27], a point cloud coding scheme is proposed. Although these point cloud compression schemes can significantly reduce the size of point cloud maps, they are more suitable for point cloud storage and transmission than for online map construction, because point cloud maps need to be compressed and decompressed frequently when this method is used for online mapping.

In some studies, multi-resolution maps are used to remove redundant information in maps. In paper [28], a multi-resolution voxel map representation method is proposed to observe the occupancy probability of different voxels, and only the voxels with higher occupancy probability are retained. This method is suitable for occupancy voxel maps. In paper [29], different resolutions are set according to the distance between the point cloud and the sensor. This method is suitable for scenarios where the sensor is fixed, but it is not easy to implement when the sensor is moving.

The variable resolution strategy proposed in this paper can consider different types of features in the map, and adopt different resolutions for different positions adaptively, which can better express the information in the environment. In addition, the map update is carried out in the local area covered by the point cloud of the current LiDAR scan, which is more convenient to implement. These features can be very helpful for some applications, such as building reconstruction and search rescue. The proposed mapping algorithm can significantly reduce the map data size during building reconstruction and search rescue while preserving adequate environmental information.

3. Method

This section introduces the entire 3D LiDAR mapping system in detail. Section 3.1 describes the overall structure of the algorithm process, Section 3.2 describes the scan-to-map based NDT (normal-distributions transform) registration algorithm, and Section 3.3 describes the variable resolution update strategy of the global map.

For clarity, some denotations used in this paper are presented as follows. The coordinate system $\left\{M\right\}$ is defined as the global coordinate system, and the point cloud in $\left\{M\right\}$ is denoted as ${}^{M}C$, where ${}^{M}C=\{{}^M{\mathbf{p}}_0,{}^M{\mathbf{p}}_1,{}^M{\mathbf{p}}_2,\cdots \}$, ${}^M{\mathbf{p}}_i\in {\mathbb{R}}^3$ is a point in point cloud ${}^{M}C$. $\left\{L_k\right\}$ is the LiDAR frame at time k, the point cloud in frame $\left\{L_k\right\}$ is denoted as ${}^{L_k}C$, where ${}^{L_k}C=\{{}^{L_k}{\mathbf{p}}_0.{}^{L_k}{\mathbf{p}}_1,{}^{L_k}{\mathbf{p}}_2,\cdots \}$, ${}^{L_k}{\mathbf{p}}_i\in {\mathbb{R}}^3$ is a point in point cloud ${}^{L_k}C$. The 3D geometry transformation from frame $\left\{L_k\right\}$ to global frame $\left\{M\right\}$ is denoted as ${}_{L_k}^M\mathbf{T}$, ${}_{L_k}^M\mathbf{T}\in SE\left(3\right)$, where $SE\left(3\right)$ is the special Euclidean group.

3.1. System Overview

The pipeline of the 3D LiDAR mapping system proposed in this paper is shown in Figure 1. The measurement error of LiDAR scanning is relatively large at a far distance, and the measured points with a small range are induced by the vehicle body itself. Therefore, a distance filter is applied to remove the aforementioned disturbing measurements. The filtered measurement is ${}^{L_k}C$. Then, the constant turn rate and velocity (CTRV) model [30] is used to correct the distortion in ${}^{L_k}C$, which results in a new point cloud ${}^{L_k}C{}^{\prime }$. Then, the reference point cloud ${}^M{C}_{ref}$ is extracted from the global map ${}^{M}C$ to overlap with ${}^{L_k}C{}^{\prime }$. A scan-to-map NDT algorithm is used to register ${}^{L_k}C{}^{\prime }$ and ${}^M{C}_{ref}$, and the transformation ${}_{L_k}^M\mathbf{T}$ from frame $\left\{L_k\right\}$ to the global frame $\left\{M\right\}$ is obtained. ${}^{L_k}C{}^{\prime }$ can be fine-tuned using the registration results, which outputs a corrected point cloud ${}^{L_k}{C}^{″}$. ${}^{L_k}{C}^{″}$ is merged to ${}^{M}C$ with an adaptive variable resolution strategy, which, as it is mentioned above, provides simplified and accurate environment models.

The global map is managed using an Octree [31] for efficient voxel queries, where the point cloud belonging to a voxel is stored in the corresponding cell. The voxel-related properties, such as the mean and covariance of the local point cloud, the curvature, and the number of the points, are stored in a hash table [32], denoted as H, for fast information extraction.

3.2. Scan-to-Map NDT Registration

Compared with the scan-to-scan registration method, the scan-to-map registration method can largely avoid cumulative registration errors [33], so the scan-to-map method is used for point cloud registration in this paper.

3.2.1. Select the Global Reference Cloud

To reduce the computational burden, the scan-to-map registration method in this paper is not to register the current point cloud ${}^{L_k}C{}^{\prime }$ directly with the global point cloud ${}^{M}C$ but to select a sub-map point cloud ${}^M{C}_{ref}$ in ${}^{M}C$ that has overlapping areas with ${}^{L_k}C{}^{\prime }$.

${}^{L_k}C{}^{\prime }$ is converted to the global frame $\left\{M\right\}$ with Equation (1), and denoted by ${}^{M}C{}^{\prime }$.

$${}^M{\mathbf{p}}^{{}^{\prime }}={}_{L_k}^M{\mathbf{R}}_{pre}·{}^{L_k}{\mathbf{p}}^{\prime }+{}_{L_k}^M{\mathbf{t}}_{pre}$$

(1)

where ${}^{L_k}{\mathbf{p}}^{\prime }\in {}^{L_k}C{}^{\prime }$, ${}^M{\mathbf{p}}^{\prime }\in {}^{M}C{}^{\prime }$. ${}^M{\mathbf{p}}^{\prime }$ are the points after ${}^{L_k}{\mathbf{p}}^{\prime }$ is transformed to frame $\left\{M\right\}$, ${}_{L_k}^M{\mathbf{R}}_{pre}\in {\mathbb{R}}^{3\times 3}$ and ${}_{L_k}^M{\mathbf{t}}_{pre}\in {\mathbb{R}}^3$ are the rotation matrix and translation vector of the predicted frame $\left\{L_k\right\}$ with respect to the global frame $\left\{M\right\}$. ${}_{L_k}^M{\mathbf{R}}_{pre}$ and ${}_{L_k}^M{\mathbf{t}}_{pre}$ are calculated from Equations (2) and (3), respectively, using the CTRV model.

$${}_{L_k}^M{\mathbf{R}}_{pre}={}_{L_{k-1}}^M\mathbf{R}·{}_{L_{k-1}}^{L_{k-2}}\mathbf{R}$$

(2)

$${}_{L_k}^M{\mathbf{t}}_{pre}={}_{L_{k-1}}^M\mathbf{R}·{}_{L_{k-1}}^{L_{k-2}}\mathbf{t}+{}_{L_{k-1}}^M\mathbf{t}$$

(3)

where ${}_{L_{k-1}}^M\mathbf{R}\in {\mathbb{R}}^{3\times 3}$ and ${}_{L_{k-1}}^M\mathbf{t}\in {\mathbb{R}}^3$ are the rotation matrix and translation vector of the frame $\left\{L_{k-1}\right\}$ with respect to the frame $\left\{M\right\}$, and ${}_{L_{k-1}}^{L_{k-2}}\mathbf{R}\in {\mathbb{R}}^{3\times 3}$ and ${}_{L_{k-1}}^{L_{k-2}}\mathbf{t}\in {\mathbb{R}}^3$ are the rotation matrix and translation vector of the frame $\left\{L_{k-1}\right\}$ with respect to the frame $\left\{L_{k-2}\right\}$.

In order to guarantee a sufficient overlap between ${}^{M}C{}^{\prime }$ and ${}^M{C}_{ref}$ in the subsequent registration process, 26 voxels $V_{adj}$ adjacent to $V_i$ in the global map point cloud ${}^{M}C$ are searched and added to the Octree of ${}^M{C}_{ref}$. Since the pose discrepancy between ${}^{M}C{}^{\prime }$ and ${}^{M}C$ is limited after applying the CTRV model, a single layer expansion with 26 adjacent voxels is enough to ensure an overlap between ${}^{M}C{}^{\prime }$ and ${}^M{C}_{ref}$. Only the mean and covariance matrix of the point clouds in voxels are useful in the NDT registration process. Therefore, we only need to record the voxel index in the Octree of ${}^M{C}_{ref}$, and the corresponding voxel properties, including the mean and covariance matrix, can be quickly found in the global hash table H.

3.2.2. The NDT Registration

The NDT is a method of compactly representing a surface. The transform maps a point cloud to a smooth surface representation, described as a set of local probability density functions (PDFs), each of which describes the shape of a section of the surface. The NDT registration is to find ${}_{L_k}^M\mathbf{T}$ that maximizes the likelihood that the points of the current scan lie on the reference scan surface [34].

The initial guess of the ${}_{L_k}^M\mathbf{T}$ is denoted as ${}_{L_k}^M{\mathbf{T}}_{pre}$, which includes ${}_{L_k}^M{\mathbf{R}}_{pre}$ and ${}_{L_k}^M{\mathbf{t}}_{pre}$. Denote the difference between the optimal solution of ${}_{L_k}^M\mathbf{T}$ and the ${}_{L_k}^M{\mathbf{T}}_{pre}$ as $\Delta \mathbf{T}$, $\Delta \mathbf{T}\in SE\left(3\right)$. The likelihood function is defined as

$$\Psi =\prod _{i=1}^n\xi \left(h(\Delta \mathbf{T},{}^M\mathbf{p}{{}^{\prime }}_i)\right)$$

(4)

or, equivalently, in its negative logistic form:

$$-log\Psi =-\sum _{i=1}^{n}log\left(\xi \left(h(\Delta \mathbf{T},{}^M\mathbf{p}{{}^{\prime }}_i)\right)\right)$$

(5)

where ${}^M\mathbf{p}{{}^{\prime }}_i\in {}^{M}C{}^{\prime }$. The definition of $h(\Delta \mathbf{T},{}^M\mathbf{p}{{}^{\prime }}_i)$ is given in Equation (6), and the definition of $\xi \left(\mathbf{p}\right)$ is given in Equation (7).

$$h(\Delta \mathbf{T},{}^M\mathbf{p}{{}^{\prime }}_i)=\Delta \mathbf{T}·{}^M\mathbf{p}{{}^{\prime }}_i$$

(6)

$$\xi \left(\mathbf{p}\right)=\frac{1}{{\left(2\pi \right)}^{3/2}{\left|{\sum }_i\right|}^{1/2}}exp\left\{-\frac{1}{2}{(\mathbf{p}-{\mu }_i)}^T{\sum }_{\mathbf{i}}^{-\mathbf{1}}(\mathbf{p}-{\mu }_i)\right\}$$

(7)

where $\xi \left(\mathbf{p}\right)$ is the probability density function corresponding to voxel $V_i$ of point cloud ${}^M{C}_{ref}$, $\mathbf{p}\in {\mathbb{R}}^3$. ${\mu }_i\in {\mathbb{R}}^3$ is the mean of points in voxel $V_i$, and ${\sum }_i\in {\mathbb{R}}^{3\times 3}$ is the covariance matrix of points in voxel $V_i$.

Points with small probability density significantly deteriorate the registration results. A mixture of a normal distribution and a uniform distribution can be used to approximate Equation (5) to eliminate the influence of outliers on the results.

$$f(\Delta \mathbf{T})=-\sum _{i=1}^{n}g\left(h(\Delta \mathbf{T},{}^M\mathbf{p}{{}^{\prime }}_i)\right)$$

(8)

where the definition of $g\left(\mathbf{p}\right)$ is given in Equation (9).

$$g\left(\mathbf{p}\right)=-d_1·exp(-\frac{d_2}{2}{(\mathbf{p}-{\mu }_{\mathbf{i}})}^T{\sum }_{\mathbf{i}}^{-\mathbf{1}}(\mathbf{p}-{\mu }_{\mathbf{i}}))$$

(9)

where $d_1$ and $d_2$ are constants. The readers are referred to the article [34] for details.

The Newton method is used to solve the following nonlinear optimization problems.

$$\Delta {\mathbf{T}}^{\ast }=\underset{\Delta \mathbf{T}}{argmin}f(\Delta \mathbf{T})$$

(10)

$$\Delta {\mathbf{T}}^{\ast }=\left[\begin{array}{cc}\Delta {\mathbf{R}}^{\ast }& \Delta {\mathbf{t}}^{\ast }\\ \mathbf{0}& \mathbf{1}\end{array}\right]$$

(11)

where $\Delta {\mathbf{R}}^{\ast }\in {\mathbb{R}}^{3\times 3}$ is the rotation matrix and $\Delta {\mathbf{t}}^{\ast }\in {\mathbb{R}}^3$ is the translation vector.

Accordingly, the rotation ${}_{L_k}^M\mathbf{R}$ and translation ${}_{L_k}^M\mathbf{t}$ of the LiDAR frame can be updated, as shown in Equations (12) and (13).

$${}_{L_k}^M\mathbf{R}=\Delta {\mathbf{R}}^{\ast }·{}_{L_k}^M{\mathbf{R}}_{pre}$$

(12)

$${}_{L_k}^M\mathbf{t}=\Delta {\mathbf{R}}^{\ast }·{}_{L_k}^M{\mathbf{t}}_{pre}+\Delta {\mathbf{t}}^{\ast }$$

(13)

3.3. Map Maintenance

With the pose of the current scan registered, the new point cloud needs to be merged into the global map, which is called global map updating. The first step is to refine the point cloud correction for ${}^{L_k}C$ with the optimized registration results. The rotation matrix ${}_{L_k}^{L_{k-1}}\mathbf{R}$ and translation vector ${}_{L_k}^{L_{k-1}}\mathbf{t}$ between the frame k and $k-1$ is calculated in Equations (14) and (15).

$${}_{L_k}^{L_{k-1}}\mathbf{R}={}_{L_{k-1}}^M{\mathbf{R}}^{-1}·{}_{L_k}^M\mathbf{R}$$

(14)

$${}_{L_k}^{L_{k-1}}\mathbf{t}={}_{L_{k-1}}^M{\mathbf{R}}^{-1}·({}_{L_k}^M\mathbf{t}-{}_{L_{k-1}}^M\mathbf{t})$$

(15)

The re-corrected scan ${}^{L_k}{C}^{″}$ is converted to the global frame $\left\{M\right\}$ with Equation (16), and denoted as ${}^M{C}^{″}$.

$${}^M\mathbf{p}{}^{″}={}_{L_k}^M\mathbf{R}·{}^{L_k}\mathbf{p}{}^{″}+{}_{L_k}^M\mathbf{t}$$

(16)

where ${}^{L_k}\mathbf{p}{}^{″}\in {}^{L_k}{C}^{″}$, ${}^M\mathbf{p}{}^{″}\in {}^M{C}^{″}$.

The essence of our variable resolution method is to maintain a geometric-adaptive point density in the map updating process. Denote the changed voxel set as V. Firstly, the curvature ${\alpha }_i$ of the point cloud in voxel $V_i\in V$ is calculated according to Equation (17).

$${\alpha }_i=\frac{{\lambda }_{0\_i}}{{\lambda }_{0\_i}+{\lambda }_{1\_i}+{\lambda }_{2\_i}}$$

(17)

where ${\lambda }_{0\_i}$, ${\lambda }_{1\_i}$ and ${\lambda }_{2\_i}$ are the eigenvalues of the covariance matrix of the point cloud in voxel $V_i$, and satisfy ${\lambda }_{0\_i}\le {\lambda }_{1\_i}\le {\lambda }_{2\_i}$.

Then, the density upper limit ${\rho }_{lim\_i}$ of the point cloud in voxel $V_i$ is calculated according to Equation (18).

$${\rho }_{lim\_i}=f_{max\_density}\left({\alpha }_i\right)=\left\{\begin{array}{cc}{\rho }_{min}& {\alpha }_i<{\alpha }_{min}\\ \eta ·{\alpha }_i& {\alpha }_{min}\le {\alpha }_i\le {\alpha }_{max}\\ {\rho }_{max}& {\alpha }_i>{\alpha }_{max}\end{array}\right.$$

(18)

where ${\rho }_{min}$ is the lower bound of ${\rho }_{lim\_i}$, ${\rho }_{max}$ is the upper bound of ${\rho }_{lim\_i}$, $\eta$ is a user-defined coefficient, ${\alpha }_{min}$ is the lower bound of ${\alpha }_i$, and ${\alpha }_{max}$ is the upper bound of ${\alpha }_i$. They can be tuned under the constraints ${\alpha }_{min}=\frac{{\rho }_{min}}{\eta }$ and ${\alpha }_{max}=\frac{{\rho }_{max}}{\eta }$. The higher the curvature of the point cloud in the voxel, the higher the upper limit of the density of the point cloud. If we need a denser map, we can set ${\rho }_{min}$, ${\rho }_{max}$, and $\eta$ larger, and vice versa.

Then, the upper limit $n_{max\_i}$ of points number in voxel $V_i$ can be calculated, as shown in Equation (19).

$$n_{max\_i}={\rho }_{lim\_i}·L^3$$

(19)

where L is the edge length of the voxel.

If the point number $n_i$ after map merging is larger than $n_{max\_i}$, it implies that redundant points exist in the $V_i$ and need to be removed. Firstly, $\gamma \ast n_{max\_i}$ points with larger probability density values are separated as the maintaining points, then $(1-\gamma )\ast n_{max\_i}$ points are randomly selected from the left part and treated as the maintaining points as well, the other points which do not belong to the maintaining points will be removed. Here, $\gamma$ is an adjustable coefficient. The probability density value of each point in voxel $V_i$ is calculated as shown in Equation (7). In the process of NDT registration, outliers will have a certain impact on the registration results, so reserving the points with high probability density will help improve the accuracy of NDT registration and then improve the accuracy of map construction.

The map maintenance process described above is concluded as in Algorithm 1. In addition, we need to update the mean ${\mu }_i$ and covariance matrix ${\sum }_i$ of the new point cloud in each voxel $V_i\in V$ in the global hash table H for further registration usage.   

Algorithm 1:Variable Resolution Strategy

4. Experiment Results and Discussion

We verified the variable resolution NDT method in three different environments (room environment, corridor environment, and outdoor environment) for performance evaluation. In addition, we carried out a long-time experiment to illustrate the advantages of our mapping algorithm in decreasing the map data size.

4.1. Equipment Introduction

In the experiment, a mobile robot equipped with a 16-line rotating LiDAR (as shown in Figure 2a) was used to conduct the environment mapping experiment. The LiDAR can rotate around the X-axis shown in Figure 2a with respect to the vehicle body to obtain a larger field of view for environmental mapping. Specific device parameters are shown in

Table 1. The parameters of the spin LiDAR equipped on the mobile robot.

Sensor Type	Line Number	FOV	Scan Frequency	Frequency around X Axis	System Accuracy
VLP-16	16	${270}^{\circ }$	10 Hz	0.2 Hz	0.006 m *

∗ System accuracy in the range of 5 m.

. During the data-capturing process, the robot navigates along the given path, and the LiDAR captures the point cloud continuously. The ground truth of the environment is captured by a high-accuracy laser scanner system that can rotate around its Z-axis, as shown in Figure 2b, whose system accuracy is 0.003 m in the range of 5 m. Specific device parameters are shown in

Table 2. The parameters of the portable laser scanner.

Sensor Type	Line Number	FOV	Scan Frequency	Rotate Range around Z Axis	System Accuracy
Hokuyo UST-30LX	1	${270}^{\circ }$	40 Hz	${180}^{\circ }$	0.003 m *

∗ System accuracy in the range of 5 m.

.

4.2. Mapping Accuracy Analysis

To verify the efficiency of the variable resolution strategy on the accuracy of mapping, we carried out experiments in the room, corridor, and outdoor environments, respectively, and collected the environment truth value to compare the point-to-point map error with the algorithm without variable resolution strategy.

4.2.1. Room Environment

The room environment mapping experiment is carried out in the scene as shown in Figure 3. The scene is a laboratory with some experimental equipment and desks placed against the wall.

The environment ground truth is shown in Figure 4(a-1,a-2), the mapping result without considering the variable resolution strategy is shown in Figure 4(b-1,b-2), and the mapping result considering the variable strategy is shown in Figure 4(c-1,c-2). The full map is on the left and the interior details are on the right. Because the LiDAR can scan through the glass to get the outdoor scene, the map is not only the indoor scene but also part of the outdoor railing. There is no point cloud on the ground in the upper right part of Figure 4(b-2,c-2), because the mobile robot moves in the lower left part of the scene when collecting the point cloud data. It can be seen from Figure 4(c-1,c-2) that the density of the point cloud in different parts of the map is different. The density of the point cloud is low in areas with small curvature (such as ceiling, wall, and floor) and high in areas with large curvature (such as the experimental equipment). This effectively reduces the number of redundant points in the map.

The CloudCompare software is used to align the point clouds shown in Figure 4(b-1,c-1) with the ground truth shown in Figure 4(a-1), and then, we randomly select some areas shown in Figure 5 to calculate the mean absolute error (MAE). The results are shown in

Table 3. The mapping error (MAE/m) comparison of the room environment.

Area Index	Without Variable Resolution	With Variable Resolution	Area Index	Without Variable Resolution	With Variable Resolution
1	0.0039	0.0031	5	0.0055	0.0052
2	0.0048	0.0041	6	0.0115	0.0112
3	0.0062	0.0061	7	0.0097	0.0093
4	0.0064	0.0063	8	0.0046	0.0042

Smaller errors are in bold.

. We selected a total of 8 areas to calculate the mapping error, among which areas 1–3 are the ceiling, areas 5 and 8 are the walls, area 4 is a piece of experimental equipment, and areas 6 and 7 are the desks. The mapping errors in these areas are smaller when the variable resolution processing strategy is used than when the strategy is not used. The accuracy of area 2 is improved most obviously, and the error is reduced by $14.58\%$.

4.2.2. Corridor Environment

The corridor environment mapping experiment is carried out in the scene as shown in Figure 6. The scene is a U-shaped corridor with walls on one side and railings on the other.

The environment ground truth is shown in Figure 7(a-1,a-2), the mapping result without considering the variable resolution strategy is shown in Figure 7(b-1,b-2), and the mapping result considering the variable strategy is shown in Figure 7(c-1,c-2). It can be seen that after considering the variable resolution strategy for the map, the density of the point cloud in the areas with low curvature (such as ceiling and floor) in the map is low, and the density of the point cloud in the areas with high curvature (such as railing) is high.

We selected a total of 16 areas (as shown in Figure 8) to calculate the mapping error (as shown in

Table 4. The mapping error (MAE/m) comparison of the corridor environment.

Area Index	Without Variable Resolution	With Variable Resolution	Area Index	Without Variable Resolution	With Variable Resolution
1	0.0743	0.0694	9	0.0787	0.0627
2	0.0304	0.0284	10	0.0657	0.0435
3	0.0418	0.0394	11	0.0243	0.0217
4	0.0569	0.0503	12	0.0382	0.0283
5	0.0977	0.0973	13	0.0571	0.0560
6	0.0207	0.0188	14	0.0729	0.0316
7	0.0186	0.0171	15	0.0176	0.0167
8	0.0682	0.0573	16	0.0623	0.0444

Smaller errors are in bold.

), among which areas 1–5 are the ceiling, areas 6 and 7 are the floor, areas 8–13 are the railings, and areas 14–16 are the walls. When considering the variable resolution strategy, the mapping errors in these regions are smaller than those without considering the strategy. The accuracy of area 14 is improved most obviously, and the error is reduced by $56.65\%$.

4.2.3. Outdoor Environment

The outdoor environment mapping experiment is carried out in the scene as shown in Figure 9. This scene is a section of the road on the campus. Several cars are parked on the road, many tall trees on the right side, a flower bed on the left, and some short trees on the flower bed.

The environment ground truth is shown in Figure 10(a-1,a-2), the mapping result without considering the variable resolution strategy is shown in Figure 10(b-1,b-2), and the mapping result considering the variable strategy is shown in Figure 10(c-1,c-2). The ground section in Figure 10(a-1,a-2) has several circular vacancies evenly spaced. The reason is that the portable laser scanner cannot scan the area under the device while collecting environmental truth values, and the scanners in adjacent positions are far apart, so they cannot cover each other’s scanning vacancies. If adjacent scannings are captured in a small distance, as in the indoor scene, the whole process would be time-consuming and inefficient. Figure 10(c-2) shows the effect of the variable resolution strategy. In the map, the part with small curvature (such as ground) has a small density of point clouds, while the part with large curvature (such as trees and cars) has a large density of point clouds.

We selected 30 areas (as shown in Figure 11) to calculate the mapping error (as shown in

Table 5. The mapping error (MAE/m) comparison of the outdoor environment.

Area Index	Without Variable Resolution	With Variable Resolution	Area Index	Without Variable Resolution	With Variable Resolution
1	0.0545	0.0417	16	0.0327	0.0327
2	0.0522	0.0416	17	0.0358	0.0356
3	0.0579	0.0517	18	0.0277	0.0287
4	0.0400	0.0210	19	0.0299	0.0290
5	0.0218	0.0194	20	0.0490	0.0485
6	0.0288	0.0254	21	0.0546	0.0541
7	0.0196	0.0188	22	0.0640	0.0630
8	0.0284	0.0261	23	0.0737	0.0735
9	0.0298	0.0269	24	0.0437	0.0423
10	0.0344	0.0340	25	0.0466	0.0453
11	0.0241	0.0234	26	0.0371	0.0359
12	0.0390	0.0382	27	0.0474	0.0456
13	0.0417	0.0374	28	0.0536	0.0516
14	0.0324	0.0294	29	0.0642	0.0621
15	0.0245	0.0218	30	0.0515	0.0499

Smaller errors are in bold.

), among which areas 1–5 are the ground, areas 6–16 are the flower bed, areas 17–19 are the cars, areas 20–23 are the tall trees, and areas 24–30 are the short trees. The mapping errors of 28 regions with the variable resolution strategy are smaller than those without the strategy. The accuracy of area 4 is improved most obviously, and the error is reduced by $47.50\%$. In the above three scenes, the most apparent error improvement is in the plane region, that is, the region with small curvature. In the variable resolution map, more abnormal points will be removed from the area with small curvature, so the area with small curvature is more likely to have a significant error improvement.

4.3. Map Size Comparison

In the three experiments in Section 4.2, the comparison of the map size is shown in

Table 6. The map size comparison.

Method	Map Size (Point Number)
	Room	Corridor	Outdoor
without variable resolution	6,173,086	10,355,562	10,986,566
with variable resolution	3,767,088	6,490,328	8,774,072

Smaller numbers are in bold.

. For the room environment, the size of the global map with variable resolution strategy is reduced to $61.02\%$ of that without variable resolution strategy. For the corridor environment, the map size is reduced to $62.67\%$. For the outdoor environment, the map size is reduced to $79.86\%$. Because there are relatively few areas with small curvature in the outdoor environment, the proportion of map size reduction is smaller than in the room and corridor environment.

4.4. Long-Time Running Experiment

In order to verify the advantages of the variable resolution map update strategy proposed in this paper to build maps in the same space for a long time, but the map does not grow infinitely over time, we conducted repeated mapping for a long time in the experimental scene of a conference room as shown in Figure 12 to cover most details in the scene. This conference room scene has several tables and chairs in the middle of the room and some chairs at the side of the room. The roof of the conference room is a strip multi-storey non-planar structure.

The environment ground truth is shown in Figure 13(a-1,a-2), the mapping result without considering the variable resolution strategy is shown in Figure 13(b-1,b-2), and the mapping result considering the variable strategy is shown in Figure 13(c-1,c-2). Since there are many flat parts in this environment (for example, the ground and walls), most of the point cloud density is small after the variable resolution strategy is applied, as shown in Figure 13(c-1,c-2). Since the ceiling in this environment is not a planar structure, the curvature is relatively high, and the point cloud density is also rather large when the variable resolution strategy is used, as shown in Figure 13(c-1).

We selected 18 areas (as shown in Figure 14) to calculate the mapping error (as shown in

Table 7. The mapping error (MAE/m) comparison of the conference room environment.

Area Index	Without Variable Resolution	With Variable Resolution	Area Index	Without Variable Resolution	With Variable Resolution
1	0.0245	0.0231	10	0.0083	0.0074
2	0.0147	0.0129	11	0.0072	0.0057
3	0.0153	0.0138	12	0.0148	0.0125
4	0.0049	0.0043	13	0.0104	0.0076
5	0.0110	0.0102	14	0.0088	0.0075
6	0.0075	0.0070	15	0.0101	0.0092
7	0.0140	0.0119	16	0.0055	0.0107
8	0.0075	0.0064	17	0.0092	0.0088
9	0.0177	0.0153	18	0.0078	0.0073

Smaller errors are in bold.

), among which areas 1–6 are the ceiling, areas 7–9 are the floor, areas 10 and 11 are the walls, areas 12–15 are the chairs, and areas 16–18 are the tables. The mapping errors of 17 regions are smaller when the variable resolution strategy is considered compared with when the strategy is not considered.

The size of the maps obtained by the two algorithms is 17,677,016 (without variable resolution strategy) and 6,438,981 (with variable resolution strategy), respectively. The size of the global map with the variable resolution strategy is reduced to $36.43\%$ of that without the variable resolution strategy. The data reduction is expected to be more significant if the experiment time is further extended.

In the mapping process of the two algorithms, the comparison of the size of the global map over time is shown in Figure 15. When we consider using the variable resolution strategy on the map during map construction, the size of the global map will almost stop growing after the scope of the map covers the experimental space rather than increase infinitely with time, as shown in the blue curve in Figure 15. When the variable resolution strategy is not considered for the map construction, the size of the global map will increase infinitely over time, as shown in the red curve in Figure 15. This advantage of the variable resolution strategy proposed in this paper is suitable for long-time mapping because the map size grows with space rather than time.

Figure 16 shows the mapping result of the three time points in the blue line in Figure 15. Figure 16a represents the mapping result at 60 s, at which most of the information in the scene was already constructed. Figure 16b shows the mapping result at 90 s, at which all the scene’s information was almost constructed. Figure 16c illustrates the mapping result at the last moment, at which the detailed information in the scene was further improved. After the 90 s, the scope of map construction covered all the space in the scene, so the map size almost stayed the same.

5. Conclusions

In order to solve the problem that the global map grows infinitely with time even in a fixed space in the process of 3D LiDAR mapping, this paper proposes a 3D LiDAR mapping algorithm based on a scan-to-map and NDT registration algorithm, which adopts a variable resolution strategy when updating the global map. The density of the point cloud in the global map adaptively changes with the curvature of different positions. Using this strategy, the size of the map can be reduced with the premise of ensuring the accuracy of the map construction so that the global map changes with space instead of growing infinitely with time.

Experiments in a room environment, corridor environment, outdoor environment, and conference room environment verify the effectiveness of the LiDAR mapping algorithm.

In the future, we may study LiDAR mapping in large-scale dynamic environments.