ERASOR++: Height Coding Plus Egocentric Ratio Based Dynamic Object Removal for Static Point Cloud Mapping

Categories of typical post-processing methods are summarized as follows. Cluster-based [6] [7] or semantic-based [8] [9] methods, which often rely on deep learning networks [15]. These methods typically require point cloud supervised labels for training, which can limit their generalization ability. Voxel-based methods [10] [11] [12] tend to be computationally expensive due to the calculation of misses or hits within each voxel. Visibility-based methods [13] [14] heavily depend on the hypothesis that dynamic objects always appear in front of static objects along a ray, rendering them invalid in scenarios where obstacles obstruct the line of sight between sensors and dynamic objects, or in open environments with few regular static objects.

As for the descriptor-based method like Egocentric RAtio of pSeudo Occupancy-based dynamic object Removal (ERASOR) [3] method, proposed by Lim et al., it leverages vertical column descriptors and could overcome the limitations of above with high running speed and a visibility-free characteristic. But it still has some problems like bad removal in some areas, especially when used in unstructured or vegetational environments as shown in the next part.

ERASOR method focuses on the problem of dynamic objects post-rejection in a set of generated maps. In that paper, Lim et al. brought egocentric descriptor into dynamic area test and proposed the major work using egocentric descriptor-based method. The main work of ERASOR includes leveraging a representation of points in vertical columns called Region-wise Pseudo Occupancy Descriptor (R-POD), a method called Scan Ratio Test (SRT) to fetch bins with dynamic objects, a static point retrieval method called Region-wise Ground Plane Fitting (R-GPF), and presenting a metrics, especially for static map building tasks. More details about the mentioned work can be found in [3].

Although ERASOR provides a promising performance on dynamic object removal, it still has some limitations, which can be summarized as follows.

• •

  Simple representation of descriptor. In the previous method, the representation of the bin area’s descriptor was simplified to only include the height difference along the Z-axis, which seems a substantial loss of valuable information.

• •

  Uncertainty and erroneous removal in areas with blocked ground. As shown in Fig.3, in some area that contains vegetation far from the scan center, or when buildings are positioned behind other objects, the limited number of scan points within the bin leads to potential inaccuracies in determining the status of the bins.

• •

  Isolated dynamic bin caused by drift or error of LiDAR data. As shown in Fig.4, certain bins may exhibit an isolated dynamic status despite the absence of any dynamic objects in that particular area, which leads to wasted computation and ineffective removal processes.

To deal with these limitations above, we modified the previous method and proposed ERASOR++, an improvement based on a height coding method.

LiDAR descriptors have been a considerably vital method to represent numerous 3D point cloud information, used in various 3D applications such as location and registration, which can be divided into two key categories, namely hand-crafted and deep-learning based methods. Focused on the traditional hand-crafted approaches, classification can be made by taking into account the extracted region, such as local, global, and hybrid descriptor [16].

Moreover, there are some global descriptors designed for specific applications, like Scan Context [17] and its expansion method [18] [19], which are primarily designed for loop detection tasks. These descriptors leverage spatial distribution or intensity information within each point block, serving as a source of inspiration for the presentation and description of LiDAR point clouds in this paper.

Raw point clouds obtained from scans and prior maps often contain an overwhelming amount of structural information, making calculations complex and challenging. Hence, it is necessary to simplify the information and adopt a novel method to describe the point cloud in the area of interest. Motivated by [3] [17] [18], we propose a novel approach by integrating the height difference and height encoding information, resulting in the Height Coding Descriptor (HCD) which offers a both concise and informative representation of the point cloud in the region of focus.

When the current scan comes, we first fetch the point cloud within a distance range, then get the submap from the prior map by fetching the point cloud in the same range. The Height Coding Descriptor is based on the fetched point cloud as presented in Fig.2.

Let $\mathcal{P}_{t}$ be the point cloud perceived at time step $t$; ${p_{t}}_{k}$ be a point in $\mathcal{P}_{t}$, i.e., ${p_{t}}_{k}\in\mathcal{P}_{t}$; and $\mathcal{V}_{t}$ be the points fetched from a specific range, also called Volume of Interest (VoI). Then, $\mathcal{V}_{t}$ can be defined and formulated as:

	$\displaystyle\mathcal{V}_{t}=\{{p_{t}}_{k}|{p_{t}}_{k}\in\mathcal{P}_{t},{\rho_{t}}_{k}<L_{max},$		(1)
	$\displaystyle H_{min}<{z_{t}}_{k}<H_{max}\},$

where ${p_{t}}_{k}=\{{x_{t}}_{k},{y_{t}}_{k},{z_{t}}_{k}\}$, and ${\rho_{t}}_{k}=\sqrt{{{x_{t}}_{k}}^{2}+{{y_{t}}_{k}}^{2}}$. Note that $L_{max}$, $H_{min}$ and $H_{max}$ are constant threshold of VoI. This formula limits the maximal radial distance boundary using $L_{max}$ and the valid height range of VoI using $H_{min}$ and $H_{max}$, which can cut down the number of points to be transformed and avoid the outlier.

As shown in Fig.2, the VoI $\mathcal{V}_{t}$ will be divided into bins according to azimuthal and radial position, known as sectors and rings. Let ${\mathcal{B}_{t}}_{ij}$ be the points of bins in $i$ ring and $j$ sector, then the points in each bin can be defined as:

	$\displaystyle{\mathcal{B}_{t}}_{ij}=\{{p_{t}}_{k}|{p_{t}}_{k}\in\mathcal{V}_{t},\frac{(i-1)L_{max}}{N_{r}}\leq{\rho_{t}}_{k}<\frac{iL_{max}}{N_{r}},$		(2)
	$\displaystyle\frac{(j-1)2\pi}{N_{\theta}}\leq{\theta_{t}}_{k}<\frac{j2\pi}{N_{\theta}}\},$

where ${\theta_{t}}_{k}=arctan2({y_{t}}_{k},{x_{t}}_{k})+\pi$, $N_{r}$ and $N_{\theta}$ respectively mean the total divisional numbers of rings and sectors.

In each bin, layer indexes based on height value are utilized to represent the height information in the descriptor. Let ${\mathcal{L}_{t}}_{ij}(\alpha)$ be the ${\mathcal{B}_{t}}_{ij}$ points in layer index $\alpha$ and designed as:

	$\displaystyle{\mathcal{L}_{t}}_{ij}(\alpha)=\{{p_{t}}_{k}|{p_{t}}_{k}\in{\mathcal{B}_{t}}_{ij},\frac{(\alpha-1)(H_{max}-H_{min})}{N_{l}}$		(3)
	$\displaystyle\leq{z_{t}}_{k}<\frac{\alpha(H_{max}-H_{min})}{N_{l}}\},$

where $N_{l}$ means the total layer number of indexes in each bin.

So the value of Height Coding Descriptor $\mathcal{D}_{t}$ as well as its difference descriptor ${{\mathcal{D}^{\mathcal{D}}}_{t}}_{ij}$ and encoding descriptor ${{\mathcal{D}^{\mathcal{E}}}_{t}}_{ij}$ of the bins in $i$ ring and $j$ sector can be defined as:

	$\displaystyle{{\mathcal{D}^{\mathcal{D}}}_{t}}_{ij}$	$\displaystyle=max\{{z_{t}}_{k}\}-min\{{z_{t}}_{k}\},\ {p_{t}}_{k}\in{\mathcal{B}_{t}}_{ij}$		(4)
	$\displaystyle{{\mathcal{D}^{\mathcal{E}}}_{t}}_{ij}$	$\displaystyle=\sum_{\alpha=1}^{N_{l}}{\mathcal{O}_{t}}_{ij}(\alpha)\cdot 2^{\alpha-1},$		(5)
	$\displaystyle{\mathcal{O}_{t}}_{ij}(\alpha)$	$\displaystyle=\begin{cases}1,&{\text{if}}\ \exists\ {p_{t}}_{k}\in{\mathcal{L}_{t}}_{ij}(\alpha)\\ {0,}&{\text{otherwise.}}\end{cases},$		(6)
	$\displaystyle\mathcal{D}_{t}$	$\displaystyle=\{{{\mathcal{D}^{\mathcal{D}}}_{t}}_{ij},{{\mathcal{D}^{\mathcal{E}}}_{t}}_{ij}\}.$		(7)

Note that the representation of ${{\mathcal{D}^{\mathcal{E}}}_{t}}_{ij}$ is defined for using bit operation to simplify computation and save storage. As shown in Fig.2, the height occupancy information in each layer can be represented completely using this encoding method with an optional number of total layers $N_{l}$, which is always set as 8 or its multiples to match computer storage unit.

Thus, the descriptor can not only show the extreme value of relative height but also represent other points’ height information with a simplified encoding method.

In previous work, Scan Ratio Test (SRT) was used to compare the current scan and its corresponding prior map descriptors, in order to found out bins with dynamic objects. SRT directly compares the values of height difference in the current scan and map scan by restricting the quotient of these two values, which can be described as:

$\displaystyle\mathcal{R}_{t}=\frac{{\mathcal{D}^{\mathcal{D}}}_{t}(Curr)}{{\mathcal{D}^{\mathcal{D}}}_{t}(Map)},$

(8)

which may cause bad removal as displayed in Fig.3. In that case, based on the proposed HCD, a novel test method called Height Stack Test (HST) is designed and utilized.

From Fig.3, we can see that the overlaps in layers above ground need to be checked, and the quantity of overlap parts demonstrates the possibility of being static. So the Height Stack Test parameter is designed as:

	$\displaystyle\mathcal{H}_{t}=[{\mathcal{D}^{\mathcal{E}}}_{t}(Curr)\wedge{\mathcal{D}^{\mathcal{E}}}_{t}(Map)]$		(9)
	$\displaystyle\wedge(\neg{Layer(Ground)}),$

where the sign $\wedge$ and $\neg$ respectively means the AND and NOT operation of bit coding, and $Layer(Ground)$ means the encoding layer with point cloud below the ground, which is described in the following part III-C. $\mathcal{H}_{t}$ represent overlap layers as well as their position above ground. Note that $\mathcal{H}_{t}$ is always turned to its number of occupied bits to be checked and limited to a low level when coding and implementing.

Thus, with restriction of the position and number of overlap layers, the height information in the middle position can be leveraged, and bad removal caused by blocked points in scan bins can be successfully avoided.

Ground Layer Test (GLT) is a proposed method that involves examining the ground layer within the HCD to determine the relative position for the HST. As the Z direction information of pose from LiDAR Odometry or SLAM is usually more uncertain than X and Y direction [20], the ground layer in each sequence is not united, which causes the necessity of searching ground layer in this system.

From the process above, each point in VoI is related to one bin and can be indexed from the ring number $i$ and sector number $j$ in ${\mathcal{B}_{t}}_{ij}$. Considering that ground points can always be fetched near the scan sensors, we use the index of each bin to iterate from near to far and get the occupied ring $i$, then check the number of points in each layer of these bins, and finally select the layer with enough most points as the ground layer in ${\mathcal{B}_{t}}_{ij}$. If more than three-quarters of bins in $i$ ring achieve the same ground layer $\gamma$, the value of $Layer(Ground)$, which encodes layers below the ground, is defined as:

$\displaystyle Layer(Ground)$

$\displaystyle=\sum_{\alpha=1}^{\gamma}2^{\gamma-1},$

(10)

where $\gamma$ is the sequence of layers with the satisfied number of ground points in the ring of nearest valid bins.

Thus, the ground layer is found and then the height descriptor can represent both absolute and relative information. Note that different from Ground Plane Fitting (R-GPF) in [3], our Ground Layer Test method only leverages a simple calculation to find the layer with ground points and runs independently before all the processes to give restrictive information for the following HST and R-GPF.

Surrounding Points Test (SPT) is a process between HST and R-GPF to cut down the wrong dynamic removal. As uncertainty always exists in LiDAR point cloud and some shelters may affect the result of bin status, there are some isolated bins in dynamic status but with no dynamic object in it, as shown in Fig.4.

This part focuses on correcting isolated bin status by searching for neighboring points, aiming to identify and rectify situations where a bin appears to be dynamic due to insufficient surrounding information. By analyzing points in the vicinity, this test ensures a more accurate assessment of the bin’s dynamic status.

For the prior map is generated by the accumulation of each scan frame, dynamic objects always have long tracks in maps. So the bins indeed with dynamic objects should be congregated as groups or lines. Under this circumstance, we proposed the Surrounding Points Test method to check each dynamic bin and make sure their status is right by searching their surrounding points. The search area can be described as:

	$\displaystyle{\mathcal{S}_{t}}_{ij}=\{{\mathcal{B}_{t}}_{pq}|(i-Range)\leq p<(i+Range),$		(11)
	$\displaystyle(j-Range)\leq q<(j+Range)\},$

where $Range$ means the point cloud searching range which can be changed according to the data sequence, and $Range=1$ is enough for most situations.

If there are dynamic bins in searching areas, the status of the target bin will remain unchanged. On the contrary, if there are no dynamic bins in searching areas, the status of the target bin will be changed to static status.

Thus, the isolated wrong removal is decreased and more static points are reserved in the output map with an increase in preservation rate.

Two parts of experiments were set to evaluate and compare the proposed algorithm, including overall system tests and ablation experiments. Overall system tests evaluated the precision and speed of the system with all methods proposed above, while ablation experiments evaluated the effect of partial methods. The initial ERASOR [3] was tested in all datasets for comparison.

According to the system framework in Fig.1, the HCD and GLT parts are preconditions of HST, and SPT can be solely added in ERASOR [3]. So the ablation experiments include the evaluation of system without HCD, GLT, and HST (in Section III-A, III-C and III-B) and system without SPT (in Section III-D). Note that they were respectively named as ++ w/o ABC and ++ w/o D in Table I.

The SemanticKITTI dataset [22] [23] was chosen for experiments, as it provides annotations for each point and labels the ground-truth dynamic objects that need to be removed, which is feasible for us to evaluate the algorithm system. Similar to the datasets used in ERASOR [3], we used the selected frames with the largest part of dynamic object in different sequences, that is, 4390-4530 frames in sequence 00, 150-250 frames in sequence 01, 860-950 frames in sequence 02, 2350-2670 frames in sequence 05, and 630-820 frames in sequence 07.

Note that the experiments were conducted in these frames above but recorded directly using the number of sequences in Table I. For a fair comparison, parameters in previous work [3] remained unchanged.

Algorithm precision was measured through the novel static map-oriented quantitative metrics in [3] called Preservation Rate (PR) and Rejection Rate (RR).

PR describes the effect of static point cloud preservation as an evaluation of bad removal, which is the ratio of the preserved static points in all static points. RR describes the effect of dynamic point cloud removal, which is the ratio of removed dynamic points in all dynamic points. Both of these metrics were calculated voxel-wise with a voxel size of 0.2. $F_{1}$ score is the combination metric of PR and RR.

The results of the previous and proposed algorithm as well as the ablation experiments are listed in Table I, in which the previous algorithm was replicated using its open source code. The visual comparison of results with main differences in red boxes is depicted in Fig.5, which provides a clear representation of the effectiveness of the algorithm. It is shown that in all sequences, our method consistently achieves extreme improvements in Preservation Rate and F1 score, with only a slight decrease in RR value. And Preservation Rate in our method can always be kept at a high level. Although in sequence 02, there is a defective PR caused by the pose uncertainty in the Z axis, the results still greatly outperform the previous work. Moreover, considering the difference in points quantity, the improvement in the preserved static points number by our algorithm far outweighs the changes in the dynamic points number.

Note that when testing in long-term datasets, there is no need to run the dynamic removal system in each frame. As an external parameter, the frame skip interval could also influence the PR and RR values of each method. We tested and compared these two methods with different interval parameters, and the results of sequences 05 and 07 are depicted in Fig.6. It shows that through the trend of PR-RR relationship, the results of ERASOR++ always display on the high-value side and win high-level precision compared with ERASOR.

That is, our proposed algorithm demonstrates superior precision compared to the previous ERASOR [3] algorithm across various sequences.

Algorithm speed was measured through the average time taken for one frame iteration during the execution of the main algorithm. We also compared the speed of our proposed method with the previous ERASOR work and other ablation systems, as shown in Table I.

On account of the main time-costing calculation part R-GPF, the count of this part decides the time cost of algorithms. As shown in Table II, although our proposed methods have some incremental calculations compared with the previous method, our algorithm demonstrates an equivalent or even faster execution time compared to the previous ERASOR algorithm, for the iteration counts of R-GPF are less than ERASOR methods. Considering that the previous ERASOR algorithm is already known to be at least ten times faster than other methods [3], our proposed algorithm, ERASOR++, also achieves a comparable level of efficiency.

This indicates that our algorithm not only improves upon the previous algorithm in terms of accuracy and precision but also maintains a comparable level of speed performance.

The results of ablation experiments are listed in Table I. From all the experiment sequences, the results of our system without each part demonstrate a decrease in Preservation Rate, which means these parts can avoid different kinds of bad removal and solve corresponding limitations.

Through these ablation experiments, the validation of each part can be proved. The Surrounding Points Test and Height Coding Descriptor as well as its corresponding test methods effectively decrease the bad removal and improve the system generalization.