A New Guaranteed Outlier Removal Method Based on Plane Constraints for Large-Scale LiDAR Point Cloud Registration

1 Introduction

Rapid advancements in 3D laser scanning and computer technology have made 3D point cloud registration a pivotal endeavor in computer vision and robotics, with its significance extending across various domains such as photogrammetry [Li et al., 2022], 3D scene reconstruction [Liu et al., 2023b], autonomous navigation [Lu et al., 2019], and virtual reality [Reyes-Aviles et al., 2022]. At its essence, point cloud registration involves the computation of a precise 6-Degree-ofFreedom (DoF) transformation T = [R|t] ∈ SE(3), serving the vital function of meticulously aligning the source point cloud P with the coordinate system of the target point cloud Q [Yan et al., 2022; Chen et al., 2022].

Accordingly, point cloud registration methods are primarily categorized into two groups: local-based and global-based methods [Ma and Wei, 2023]. Among these, the Iterative Closest Point (ICP) [Besl and McKay, 1992], a well-known local-based method with several variants [Segal et al., 2009; Bouaziz et al., 2013; Zhang et al., 2021; Koide et al., 2021], has had a profound influence on subsequent research in this field. Nevertheless, ICP and its variants employ a greedy and exhaustive nearest neighbor (NN) search during each iteration to establish point correspondences [Lim et al., 2022]. Consequently, ICP-variants heavily depend on the initial alignment and is susceptible to becoming trapped in local minima. Local registration methods are most effective when the source point cloud and the target point cloud are in close proximity, nearly overlapping [Pomerleau et al., 2013].

∗Corresponding author

In response to the limitations of ICP and its variants, global-based registration methods have gained attention for providing an initial alignment that aids subsequent local registration, guiding estimations towards the global minimum [Lim et al., 2022]. This paper mainly focuses on correspondence-based global registration methods, categorizing them into point and primitive correspondence-based approaches. Point correspondence-based methods typically involve two essential stages: establishing correspondences based on local point features and estimating transformation from identified correspondences [Li et al., 2023]. During establishing point correspondences, it is inevitable that some noise correspondences are introduced. When dealing with noisy correspondences, a commonly employed strategy is to utilize "correspondences + RANSAC (Random Sample Consensus) [Fischler and Bolles, 1981]" [Barath et al., 2020; Barath and Matas, 2018]. Although RANSAC-based methods are renowned for their simplicity and efficiency, it is important to note that their runtime exponentially escalates as the outlier rates increase. Another widely recognized strategy is "correspondences + outlier removal" [Bustos and Chin, 2017; Yang et al., 2020; Yan et al., 2022], which effectively filters out outlier correspondences. Nonetheless, methods relying on point features encounter challenges such as longer computation times, sensitivity to point density, and reduced descriptor discriminability due to different viewpoints.

Recently, there has been growing interest in the application of learning-based techniques in point cloud registration. On one hand, specific deep learning approaches have been dedicated to enhancing distinct components of the registration process, including the development of more discriminative point feature descriptors [Huang et al., 2021; Ao et al., 2021; Wang et al., 2022; Wang et al., 2023] and more efficient correspondence selection methods [Aoki et al., 2019; Bai et al., 2021; Fu et al., 2021; Xue et al., 2024]. On the other hand, [Li and Harada, 2022; Yew and Lee, 2022; Qin et al., 2022; Liu et al., 2023a] have showcased end-to-end point cloud registration methods. However, learning-based approaches require substantial training datasets and continue to face challenges related to generalization, especially when dealing with outliers and cross-dataset conditions. When handling large-scale outdoor environments, learning-based methods face limitations due to small trainable datasets and the large-scale of point cloud data, hindering their full potential in point cloud registration tasks.

Unlike point-based methods, primitive-based registration approaches leverage higher-level geometric features in scenes, such as lines [Yan et al., 2023; Liu et al., 2023c], planes [Chen et al., 2019; Li et al., 2022; Ma et al., 2023], and spheres [Zhang et al., 2022], making them better suited to handle challenges like uneven point density, occlusions, and noise in point clouds. While existing primitive-based registration methods [Chen et al., 2019; Li et al., 2022] have showcased remarkable performance on large-scale LiDAR point clouds, it is important to note that primitive descriptors often exhibit lower discriminability compared to point descriptors. This reduced discriminability can lead to a higher likelihood of introducing outlier correspondences during the establishment of primitive correspondences. Unfortunately, current primitive-based algorithms lack effective mechanisms for filtering out these outlier correspondences, which in turn restricts their precision in point cloud registration tasks.

In this paper, we present an efficient framework that decouples rotation estimation and translation estimation. For rotation estimation, our method introduces a pioneering outlier removal method based on plane primitives. For translation estimation, we propose an innovative approach for estimating relative translation using plane translation components.

In summary, the contributions of our work are as follows:

1. We present a novel point cloud registration framework grounded in plane constraints, effectively decoupling rotation estimation and translation estimation.
2. We propose a guaranteed outlier removal method for rotation estimation. This method combines rotationinvariant constraints with computational geometric consistency checks for outlier removal and globally searches for the optimal relative rotation.
3. We propose a component-wise method for translation estimation. This approach efficiently estimates each translation component through plane translation constraints and globally estimates the optimal relative translation.
4. Experimental results on multiple datasets convincingly demonstrate the robustness, effectiveness, and efficiency of our proposed framework, achieving state-of-the-art performance.

2 Methodology

First, we represent the source and target point clouds captured by a 3D LiDAR sensor as P and Q, respectively. Drawing inspiration from [Li et al., 2022], we employ a voxel-based plane extraction method to obtain planes from P and Q. Let us define P = {P1, P2, · · · , PM} and Q = {Q1,Q2, · · · ,QN}, where each plane Pi = (xi, pi)(1 ≤ i ≤ M) from P and each plane Qj = (yj, qj)(1 ≤ j ≤ N) from Q respectively contain normals xi and yj, along with points pi and qj on their respective planes. Next, let (Pk,Qk) ∈ A(1 ≤ k ≤ M ∗ N) represent candidate plane correspondences, where Pk ∈ P and Qk ∈ Q. While it is inevitable that the candidate set A contains some noisy correspondences, our goal is to minimize the impact of these outliers and accurately estimate the relative rotation R ∈ SO(3) and relative translation t ∈ R3 from candidate set A.

As shown in Figure 1, after extracting planes from the source cloud P and target cloud Q, along with the establishment of plane correspondences A, we propose a robust decoupled point cloud registration method that effectively decomposes rotation estimation and translation estimation.

2.1 Outlier Removal for Rotation Estimation

In this section, we present a guaranteed outlier removal method for estimating relative rotation from plane normals. For each candidate normal correspondence (xk, yk) from (Pk,Qk), we initially apply coarse filtering to remove outliers based on rotation-invariant constraints. Subsequently, we propose a strategy based on computational geometric consistency checks to finely prune outliers while calculating rotation bounds for each inlier (xi, yi) from (Pi,Qi) (i ̸= k) to the candidate (xk, yk). We then utilize an interval stabbing approach to estimate the rotation Rk for the candidate (xk, yk). Finally, we employ a global search to estimate the optimal relative rotation R∗ fromRk for each (xk, yk).

Rotation-Invariant Constraints

Given two pairs of non-parallel normal correspondences, the candidate (xk, yk) from (Pk, Qk) and an inlier (xi, yi) from (Pi,Qi), a rotation can be determined [Chen et al., 2019]. In an ideal situation, discounting noise, and assuming that rotation ˆR= R(ˆ θ, ˆa) precisely aligns xk to yk and xi to yi:

where ˆ θ represents the rotation angle, while ˆa denotes the rotation axis. As stated in (1), it is derived that:

Assuming that the rotation ˜R = R(˜ θ, ˜a) aligns yi to yk precisely, it is inferred that:

According to (3), one can infer that:

In the usual case involving noise and outliers, derived from (2) and (4), the rotation-invariant constraints are defined as:

where ξ represents the angular threshold for the constraints. Moreover, to avoid parallelism between the candidate (xk, yk) and an inlier (xi, yi), we set ∠(yk, yi) ∈ [ 1 12 π, 11 12 π] and ∠(xk, xi) ∈ [ 1 12 π, 11 12 π] [Li et al., 2022].

Rotation Bounds

Rotation-invariant constraints serve as a preliminary method to filter out outlier correspondences. Further, a robust and refined outlier removal method is presented in this section, enhancing the accuracy in pruning outlier correspondences.

For each candidate normal correspondence (xk, yk) and an inlier (xi, yi), our goal is to estimate Rk which maps xk to yk and brings xi within angular distance ε from yi:

As illustrated in Figure 2, the rotation matrix Rk is decomposed into two rotations:

Bk is a rotation matrix that maps xk to yk:

where [·]× is a skew-symmetric cross-product matrix.

Aθ,yk in (7) is a rotation about the axis yk by an angle θ, maintaining yk and Bkxk unchanged. It is denoted as:

where exp(·) is the exponential map [Yan et al., 2022]. Combining (8) and (9), the constraint (6) is rewritten as:

As illustrated in Figure 2, rotating by Aθ,yk while maintaining (10) causes Bkxi to be situated within a spherical region with an angular radius ε centered at yi, denoted as:

The trace of Bkxi in Figure 2 rotating around the axis yk forms a circle, which is denoted as:

As depicted in Figure 2, if circ(Bkxi, yk) and Sε(yi) do not intersect, (xi, yi) cannot be aligned by any rotation Rk and it is removed safely. If they intersect, (xi, yi) is an inlier for (xk, yk), and points pa and pb are the bounds of rotation Aθ,yk . Further, The bounds are denoted as an angle interval θ ∈ [∠(Bkxi, pa), ∠(Bkxi, pb)] = [θa i , θb i ].

As illustrated in Figure 3, the computation for the rotation bounds [θa i , θb i ] is demonstrated by projecting the sphere from Figure 2 onto a 2D plane. To simplify, we align yk with the z-axis. Let Bkxi = [x0, y0, z0]⊤. The center of circ(Bkxi, yk) is denoted as x′i = [0, 0, z0]⊤, with a radius rx′i = √x20 + y2 0 . The equation for circ(Bkxi, yk) is:

where p = [x, y, z0]⊤ and z0 is a constant value from Bkxi. Another circle, defining the boundary of Sε(yi) with an angular radius ε centered at yi, is denoted as:

The center of circ(ε, yi) is y′i = cos(ε)yi = [x1, y1, z1]⊤, and its radius ry′i = sin(ε). Point p = [x, y, z]⊤ on the circ(ε, yi) satisfies:

To obtain the intersections between the two circles circ(Bkxi, yk) and circ(ε, yi), combining (13) and (15), and letting p = [x, y, z0]⊤, the equation for line papb is:

Further, (13) and (16) are combined to obtain pa and pb:

If ∆ < 0, it indicates that the correspondence (xi, yi) is an outlier and it is removed safely. Otherwise, let pa = [ -2uv-√∆ 2(1+u2 ) , -2uv-√∆ 2(1+u2 ) u + v, z0]⊤, and pb = [ -2uv+√∆ 2(1+u2 ) , -2uv+√∆ 2(1+u2 ) u + v, z0]⊤. The rotation bounds are obtained as θi ∈ [θa i , θb i ](θa i ≤ θb i ), where θa i = ∠(Bkxi, pa), and θb i = ∠(Bkxi, pb).

Interval Stabbing

After pruning outlier correspondences for candidate (xk, yk), a set of angle intervals Θk = {[θa i , θb i ]} is obtained. An interval stabbing technique [Bustos and Chin, 2017; Cai et al., 2019] is introduced to determine an angle θk which is encompassed by maximum number of angle intervals within Θk:

where I(·) is an indicator function, yielding 1 when the input is selected and 0 otherwise. Finally, Aθk,yk Bk = Rk ∈ Rc, where Rc denotes the rotation set.

Global Rotation Search

After obtaining Rk ∈ Rc for each (xk, yk), we globally search for the optimal relative rotation R∗ :

Moreover, with the optimal R∗ , for each correspondence (Pk,Qk), if ∠(yk,R∗xk) > 2ε, (Pk,Qk) is treated as an outlier and subsequently removed from A. The experiment shows that relaxing the threshold from ε to 2ε, without affecting R∗ , enhances the precision of relative translations.

2.2 Component-Wise for Translation Estimation

Finally, a component-wise method based on plane constraints, named as CWTE, has been proposed for estimating relative translations. We define that parallel planes represent an identical translation component, and each translation component determines 1-DoF. In theory, at least three non-parallel translation components are required to determine a 3-DoF relative translation [Li et al., 2022].

We parameterize planes using the Hessian form π = [n⊤, d]⊤, where n = [nx, ny, nz]⊤ is the unit vector representing the plane orientation and d is the distance of the plane from the origin [Hartley and Zisserman, 2003]. A point p on the plane π gives d = -n⊤p. Distance between two parallel planes πi = [n⊤, di]⊤ and πj = [n⊤, dj ]⊤ is defined as:

With the optimal R∗ , for each correspondence (Pk, Qk) ∈ A, where Pk ∈ P and Qk ∈ Q, the normal xk from Pk is aligned to the normal yk fromQk by R∗ , denoted as follows:

Translation Components Computation

First, based on the parallel characteristics among planes in the source point cloud S, the plane correspondences in A′ are divided into L component sets, each component set is defined as C l (1 ≤ l ≤ L). For the l-th component set C l , an optimal 1-DoF translation component ˆtl will be estimated based on plane translation constraints.

Assuming that there are two parallel inlier plane correspondences (P ′k ,Qk), (P ′ i ,Qi) ∈ A′ , and Qk ∥ Qi. As shown in Figure 5(a), in the ideal case, the plane translation operator [0⊤, ˆ d]⊤ precisely aligns plane P ′k to plane Qk, and plane P ′ i to plane Qi, denoted as:

As indicated by (23), it is derived that:

Moreover, assuming that plane Qi is aligned to plane Qk by the plane translation operator [0⊤, ˜ d]⊤, denoted as:

According to (25), one can infer that:

In the usual case involving noise and outliers, derived from (24) and (26), the plane translation constraints are defined as:

where η represents the distance threshold for the constraints. According to the plane translation constraints in (27), we estimate the l-th translation component ˆtl , defined as follows:

where Project(π, p) denotes a function that projects the point p onto the plane π.

Global Translation Search

In accordance with (28), we have obtained the translation component set Tc = {ˆtl }Ll=1 , and each component ˆtl ∈ Tc is associated with a plane correspondence, denoted as (P ′ l ,Ql). As depicted in Figure 5(b), three different translation components ˆtu ,ˆtv ,ˆtw ∈ Tc determine a 3-DoF relative translation ˆt ∈ R3 , denoted as follows:

and we compute each candidate relative translation ˆt ∈ Tr from Tc. A global strategy is utilized to determine the optimal relative translation t∗ ∈ R3 by the following function:

3 Experimental Results

3.1 Experimental Environment

Our proposed method has been evaluated on three widely recognized LiDAR point cloud datasets, with each point cloud comprising approximately 5 to 30 million points: ETH3D [Schops et al., 2017], ETH-TLS [Theiler et al., 2015], and WHU-TLS [Dong et al., 2020]. These experiments were conducted on a PC equipped with an Intel Core i7-10700 CPU @ 2.90 GHz and 16 GB RAM.

Baselines. We conducted a comprehensive performance comparison between our method and several well-known algorithms, including FGR [Zhou et al., 2016], VGICP [Koide et al., 2021], PLADE [Chen et al., 2019], Teaser++ [Yang et al., 2020], VPFBR [Li et al., 2022], PointDSC [Bai et al., 2021], and RoReg [Wang et al., 2023]. PLADE and VPFBR represent two prominent plane-based methods, while Teaser++ and FGR are point correspondence-based outlier removal methods. PointDSC and RoReg are popular learningbased methods, and VGICP is a recently proposed variant of ICP. Notably, the learning-based methods PointDSC and RoReg were executed on a server equipped with an Intel(R) Xeon(R) 4210R CPU @ 2.40 GHz, 32GB RAM, and a GeForce RTX3090 GPU with 24GB RAM.

Metrics. We have focused on two key evaluation metrics to assess the performance of our proposed method: computational efficiency and registration accuracy. Computational efficiency is gauged by measuring runtime, while registration accuracy is quantified through rotation error and translation error [Yang et al., 2020], as expressed by the equation:

where R ∈ SO(3) and t ∈ R3 denote rotation matrix and translation vector, respectively.

Parameter Configuration. When extracting planes, we choose the top 20 planes (N = M = 20) from each point cloud, based on [Li et al., 2022]. In rotation estimation, the inlier threshold ε defaults to 1.5◦ , while the rotation-invariant constraints threshold ξ is set to 8◦ . In translation estimation, the plane translation constraints threshold η is 0.1m.

• (a) Rotation performance.
• (b) Translation performance.

Figure 6: The statistical results of the experimental rotation errors and translation errors from three datasets.

3.2 Performance Comparison

Registration Accuracy. We first evaluated the registration accuracy of our proposed algorithm. Tables 1 and 2 list the rotation and translation errors for each experiment of all algorithms on the three datasets, and Figure 6 provides a visual representation of the corresponding statistical outcomes. Evidently, our method outperformed others, exhibiting the lowest rotation and translation errors, a perspective further substantiated by the results in Table 3. VPFBR, a method based on plane pairs, ranks second in performance after our method. Conversely, another plane-based method, PLADE, despite performing well in several scenes, reveals instability, with rotation errors exceeding 5◦ in over 53% of cases. Teaser++ is a classic point correspondence-based outlier removal method, exhibiting commendable performance with rotation errors below 5◦ and translations under 1m across eight different scenes from three datasets. However, its performance significantly varies in four other scenes, resulting in rotation errors of approximately 170◦ and translation errors exceeding 20m. This notable inconsistency in results highlights the lack of robustness in the Teaser++ algorithm and underscores its varying reliability across different scenes. Unfortunately, as depicted in Figure 4, both learning-based methods, PointDSC and RoReg, exhibited poor performance in both rotation and translation. In rotation estimation, our method ensures accuracy by utilizing an outlier removal scheme, which ensures that the rotation error does not exceed ε (default ε = 1.5◦ ), as shown in Figure 2. In translation estimation, plane translation constraints facilitate efficient computation of each component, while a global search guarantees robust and accurate estimation.

• (a) Rotation performance.
• (b) Translation performance.

Figure 7: Ablation study for inter threshold ε.

Computational Efficiency. In terms of computational efficiency, as depicted in Table 3, the two learning-based methods had the shortest runtime when training time was excluded. Our proposed method not only excelled in registration accuracy but also closely followed the learning-based methods in runtime performance. Significantly, our algorithm's runtime is primarily dedicated to preprocessing. As indicated in Table 4, without parallel optimization, the average runtime for the sum of rotation and translation estimations in our method stays below 0.3s. Table 5 shows experimental results comparing our algorithm with Teaser++. Teaser++ demands a large input payload of about 1.8 × 106 bytes, including 5184 FPFH features from the source point cloud and 8584 FPFH features from the target point cloud, far exceeding our algorithm's input requirements. Conversely, our method efficiently obtains higher-precision rotation and translation from a mere 960 bytes of plane features, highlighting our algorithm's remarkable efficiency and accuracy.

3.3 Ablation Study

The inlier threshold ε is a critical parameter in our algorithm, as depicted in Figure 2, directly influencing both rotation and translation accuracy. We evaluated the impact of ε on our algorithm by selecting two scenes from each of the three datasets and varying ε = [0.5◦ , 1.0◦ , 1.5◦ , 2.0◦ , 5.0◦ ]. Figure 7 illustrates that with an inlier threshold ε between 1.0◦ and 2.0◦ , rotation error remains stable, generally below 1.0◦ , and translation error stays under 1m. Notably, extremely small ε (≤ 0.5◦ ) may overly suppress outliers, potentially reducing rotation and translation accuracy. Conversely, excessively large ε (≥ 5.0◦ ) may compromise outlier removal, adversely affecting accuracy, too.

4 Conclusion

In this paper, we introduce an innovative plane constraintsbased framework for decoupled point cloud registration in large-scale LiDAR point cloud scenes. For rotation estimation, we present a robust outlier removal scheme based on rotation-invariant constraints and computational geometric consistency checks to estimate relative rotations accurately. For translation estimation, we propose a componentwise method based on plane translation constraints to obtain high-precision relative translations. In the future, we will explore the application of our method in multi-view point cloud registration tasks.

Acknowledgments

This work was supported by the National Science Foundation of China under Grant 61771146.