Pairwise coarse registration of point clouds by traversing voxel-based 2-plane bases

ABSTRACT

To obtain the complete coverage when scanning a large-scale scene using a laser scanner, registration of pairwise point clouds by the corresponding geometric features is usually necessary under the condition of insufficient direct geographic reference information of GNSS sensors. In this paper, we propose a target-less coarse registration method for pairwise point clouds using the corresponding voxel-based planes, which are identified by a pair of conjugate 2-plane bases. The voxel-based planes are firstly extracted from the entire point cloud with 3D cubic grids decomposed by the octree-based voxelization; Then, the 2-plane bases in each point cloud are constructed, and by comparing the dihedral angles of two 2-plane bases that come from the source and target point clouds, respectively, the conjugate 2-plane base pairs are generated one by one; Next, a set of plane correspondences is identified by a pair of conjugate 2-plane bases, and its corresponding largest consistency planes (LCP) set is calculated; Finally, a series of plane correspondence sets are obtained using the generated pairs of conjugate 2-plane bases, and the one with the highest LCP is used to compute the transformation matrix between the pairwise point clouds. Experimental results revealed that our proposed pairwise coarse registration method can be effective for aligning point clouds acquired from indoor and outdoor scenes, with rotation errors less than 0.4 degrees, translation errors less than 0.4 m, root mean square distance (RMSD) less than 0.42 m, and successful registration rate (SRR) about 98%. Furthermore, our proposed method was more efficient than the point- and line-based methods under the same hardware and software configuration conditions.

KEYWORDS:

• pairwise coarse registration
• voxel-based planes
• 2-plane bases
• largest consistency planes

1. Introduction

The light detection and ranging (LiDAR) technology has been developed rapidly over the past few decades, which can obtain millions of 3D points in a few seconds, and these data have become widely used in an extensive variety of applications, such as urban planning (Gehrung et al. 2016; Li et al. 2017), scene reconstruction (Choi, Zhou, and Koltun 2015; Dong et al. 2018a), digitization management of cultural heritage (Remondino 2011; Yang and Zang 2014; Sánchez-Aparicio et al. 2018), deformation monitoring (Awrangjeb, Gilani, and Siddiqui 2018; Chen, Yu, and Wu 2018; Kusari et al. 2019; Zang et al. 2019), and forest monitoring (Kelbe et al. 2016; Polewski et al. 2018). Due to the limited view of a laser scanner, a single scan is not sufficient to describe the entire target area. Therefore, multiple scans need to be connected to get the whole information of the measured object; That is, all individual point clouds should be aligned to a common coordinate system (Dong et al. 2018b), which is an essential technology in 3D data processing (Wang et al. 2021). This process includes two steps: coarse and fine registration. Fine registration usually adopts the well-known iterative closet point (ICP) (Besl and McKay 1992) algorithm or its variants (Li and Song 2015; Tao et al. 2018; Shi, Liu, and Han 2020). However, the ICP-based methods usually require an appropriate initial transformation matrix and their iterative process is time-consuming. Without a coarse registration as a priori, this kind of methods tends to converge to a local minimum or does not converge at all. Therefore, coarse registration is of vital importance, which is just to calculate the transformation matrix between two scans only using the geomatic information that is obtained from the point cloud directly but without any additional complements. In this paper, our goal is to perform the coarse registration based on voxel-based plane congruent sets. Many researchers have developed a large number of methods to register pairwise point clouds roughly. These methods compute the transformation matrix with six degrees of freedom (DOF) (i.e. three translations and three rotations) by extracting different geometric features, including point, line, and plane.

The point-based approaches are the most common and many algorithms have been proposed. Some methods first extract the keypoints from pairwise point clouds to be matched, and then calculate the corresponding local shape descriptors. By comparing the similarity of these descriptors, the correspondences between keypoints from different point clouds can be established, which are used to compute the transformation matrix. These methods are called descriptor-based methods. Rusu, Blodow, and Beetz (2009) proposed the fast point feature histogram (FPFH) descriptor, and the possible correct corresponding points were established using the sample consensus method by comparing their FPFHs, then the transformation matrix was obtained. Fu et al. (2021a) designed a novel descriptor called multiscale eigenvalue statistic-based descriptor (MEVS), which was then used to perform the coarse registration for mobile laser scanning (MLS) point clouds. Based on the covariance matrix of the keypoints, Fu et al. (2021b) constructed the feature descriptor vector (FDV), and the pairwise terrestrial laser scanner (TLS) point clouds were aligned. Besides the above-mentioned descriptors, a lot of other descriptors have also been reported in the literatures, including quintuple local coordinate image (QLCI) (Tao et al. 2020a), binary shape context (BSC) (Dong et al. 2018b), rotational contour signature (RCS) (Yang, Xiao, and Cao 2019), local voxelized structure (LoVS) (Quan et al. 2108), and so on. The descriptor-based methods can well align pairwise point clouds with large overlap, but they are easy to fail when the overlap is small. There are also some other methods that are different from descriptor-based methods. For example, the four-point congruent sets (4PCS) (Aiger, Mitra, and Cohen-Or 2008) also belong to the point-based methods, which utilized the rule of intersection ratios between four points to establish correspondences and calculate the transformation matrix together with the registration score. By repeating establishing correspondences and calculating transformation as well as its score, the transformation matrix with the highest score was regarded as the final registration result. In order to improve the registration accuracy and efficiency, some improvements were made for 4PCS, e.g. Super4PCS (Mellado, Aiger, and Mitra 2014), keypoint-based 4PCS (K-4PCS) (Theiler, Wegner, and Schindler 2014) and semantic keypoint-based 4PCS (SK-4PCS) (Ge 2017). The 4PCS-based algorithms are robust to noise, outliers, and small overlap between pairwise point clouds to some extent. But their registration accuracy and efficiency need to be further improved. There are also some other point-based algorithms besides the before-mentioned methods, such as coherent point drift (CPD) algorithm (Myronenko and Song 2010) and its variants (Lu et al. 2016; Zang et al. 2020). This kind of methods solves the transformation matrix between pairwise point clouds using the method of probability density estimation, but which cannot well handle point clouds with large volume and are easily affected by noises and outliers. Moreover, as a promising method, deep learning-based methods (Aoki et al. 2019; Zhang et al. 2019; Pujol-Miró, Casas, and Ruiz-Hidalgo 2019) can learn advanced information from mass data to well register pairwise point clouds with large overlap and small-scale, but the result of aligning a pair of point clouds with small overlap and large-scale is poor.

In a word, the point-based methods have good versatility and can obtain enough good registration results for large overlapped point clouds, but they possess high computation complexity and low robustness to noises, occlusions and limited overlap.

Compared with point features, lines are higher level geometric features with better robustness in surveying and remote sensing. Many scholars have studied the registration method based on line features. Al-Durgham et al. (2013, 2014) first intersected the planes separated from the origin point cloud to extract 3D line features; Then, the included angle and distance between all pairs of lines were calculated and the candidate corresponding lines were selected; Finally, the transformation matrix was calculated through the weight correction process (Renaudin, Habib, and Kersting 2011). The above two methods (Al-Durgham, Habib, and Kwak 2013; Al-Durgham and Habib 2014) use non-conjugate point pairs of the conjugate 3D lines to compute the transformation parameters, although adjusted by the weight process, the accuracy is still poor. Yang et al. (2015) extracted and matched 2D building contours to register different types of urban point clouds. Their method provides a robust and reliable solution for registering TLS and airborne laser scanning (ALS) data, but it is only suitable for specific datasets with many buildings. Ge and Hu (2020), Li et al. (2021a) and Tao et al. (2020b) first projected point clouds along the vertical (z-axis) direction to extract 2D line features; Then, constructed 2D line structural unit pairs for horizontal rough registration; Finally, calculated the vertical alignment and obtained the 3D transformation matrix. These three methods (Ge and Hu 2020; Li et al. 2021a; Tao et al. 2020b) can only well register pairwise point clouds that there are only rotation differences along the z-axis. Xu, Xu, and Yang (2022) proposed a coarse registration method for urban point clouds based on 2-lines congruent sets (2LCS). First, the vertical and horizontal direction lines were extracted to construct the 2-lines bases; Then, according to the geometric character, the congruent 2-lines base pairs were identified and the 2LCS was formed; Finally, transformation parameters were estimated using the 2LCS. While determining conjugate 2-lines bases, the length and midpoint of the line were adopted. Due to shielding and difference of scanning views, the length of a line may be inconsistent in different scanning, which is not invariant character between pairwise point clouds.

Generally speaking, line-based methods can well register pairwise point clouds to some extent, and are more robust and efficient than point-based methods on aligning point clouds with large-scale and small overlap. However, 2Dline-methods can only register point clouds obtained by a strictly levelled TLS, because which hypothesizes that there only exist rotation differences along the z-axis; 3Dline-based methods lack the support of theory, because the points on the corresponding line pairs cannot be guaranteed a pair of conjugate points.

Regarding plane-based methods, the four parameters (a unit normal vector and a constant term denoting the distance from the origin to the plane) of a plane are unique, even obtained from different scans and partly. Therefore, this kind of methods have the theoretical support while solving the transformation matrix by the parameters of a pair of corresponding planes. Moreover, plane features are usually less than point or line features in a point cloud, which indicates that plane-based methods are more time-saving. Thus, the plane-based registration methods are worthy to be studied. Li et al. (2021b; 2022) applied the relationship between point and plane (Li et al. 2021b), line and plane (Li et al. 2022) as the constraint, respectively, and realized the accurate registration between pairwise point clouds using the duality four elements method. However, the two methods rely on human-computer interaction to select the corresponding primitives, which is difficult to be widely applied in the field of remote sensing due to its poor degree of automation. Xiao et al. (2013) first segmented the planes from the original point cloud using the cache octree region growth method; Then, two planes with similar areas from source and target point clouds were selected as a pair of potential corresponding planes to calculate the transformation matrix. Liang et al. (2016) first grouped the extracted planes according to the angle between their unit normal vectors; Then, the corresponding plane pairs were identified from the pairwise point clouds by comparing their dihedral angles and plane areas; Finally, the plane correspondences were used to calculate the transformation matrix. Favre et al. (2021) introduced a new metric based on plane characteristics (e.g. the distance between the projections of the origin on source and target planes, the distance between the centroids of source and target planes, the area ratio between the source and target planes, and the dot product of the normal vectors of the source and target planes) to find the best plane correspondences; Then, the optimal transformation was estimated using the plane correspondences through a two-step minimization approach: closed-form optimization and Gauss-Newton minimization. The before-mentioned three methods (Xiao et al. 2013; Liang et al. 2016; Favre et al. 2021) all use the plane area as a constraint for correspondences identifying, which requires to keep the completeness of a plane while scanning. However, due to the scanning shade, it is less possible to obtain one plane completely and with the same areas in different scans. Xu et al. (2019) first extracted planar patches from voxelized point clouds by merging plane voxels; Then, the ratios of plane normal vectors were used to find the candidate corresponding four planar patches from the pairwise point clouds; Finally, the voxel-based four-plane congruent sets (V4PCS) method was proposed to register pairwise point clouds. This method just uses plane parameters to identify the correspondences and realize registration, which can obtain high alignment accuracy. But the process of generating V4PCS bases in each point cloud and constructing conjugate V4PCS bases from different point clouds are complicated.

To overcome the high complexity of the pairwise point clouds registration, as well as guarantee its theoretical foundation and efficiency, in this paper, a pairwise coarse registration method by traversing voxel-based 2-plane bases was proposed, which was inspired by Xu et al. (2019). First, the point cloud was voxelized by the octree structure and these with planar properties were preserved, but without merging like what Xu et al. (2019) did; Second, the corresponding plane pairs were identified by a pair of conjugate 2-plane bases that have the similar dihedral angle; Finally, a series of transformation matrixes were computed by the corresponding plane pairs using the traversing strategy, and the optimal one was selected to register the pairwise point clouds. The main contributions of our work are summarized as follows:

1. A voxel-based plane extraction method is proposed. The points in the voxel with the planar properties are used to fit the parameters of the voxel-based plane, which will form the basic plane primitives for the proposed registration method.

2. A 2-plane bases construction and conjugate 2-plane base pairs generation algorithms are developed only by the dihedral angle. And a closed-form solution of transformation estimation using all corresponding plane pairs is adopted to optimize the coarse registration matrix.

3. A largest consistency planes (LCP) criterion is proposed to evaluate the transformation matrixes and used to select the optimal one.

The remainder of this paper is organized as follows. In Section 2, the detailed explanation of the proposed registration method is provided; Section 3 presents the experiments and an explicit discussion and analysis of the derived results; And the paper is concluded in Section 4, as well as the plan for future research is exhibited.

2. Proposed registration method

Our proposed coarse registration method consists of three core steps: (1) voxel-based planes construction, (2) conjugate 2-plane base pairs generation, and (3) registration using corresponding planes.

Step1: The point cloud of two scans is voxelized into the 3D grid structure, respectively. For each voxel, the smoothness of points within the voxel is estimated. Only those voxels that satisfy the criteria of curvature are deemed as planar voxels, and their plane parameters are calculated to form the voxel-based planes.

Step2: A 2-plane base is formed by two planes that satisfy the dihedral angle constraint in each point cloud. The conjugate 2-plane base pairs that have the similar dihedral angle from the pairwise point clouds are used to calculate the rotation matrix $R$, by which the initial set of corresponding planes $C_{ini}$ is obtained.

Step3: The initial transformation matrix $T_{ini}$ is estimated using $C_{\mathrm{i}\mathrm{n}\mathrm{i}}$ and its corresponding LCP is calculated; Besides, $C_{\mathrm{i}\mathrm{n}\mathrm{i}}$ is updated by $T_{ini}$, obtaining the optimal set of corresponding planes $C_{opt}$.

With a pair of conjugate 2-plane bases, an initial transformation matrix $T_{ini}$ and its corresponding LCP can be obtained. All 2-plane bases in the source point cloud are traversed to find their conjugate 2-plane bases from the target point cloud, then a series of initial transformation matrixes can be calculated as well as their LCPs. The matrix with the highest LCP is optimized using its corresponding optimal corresponding planes set $C_{opt}$, and the optimal matrix $T_{opt}$ is regarded as the final coarse registration result. In Figure 1, the processing workflow is sketched.

Figure 1. Workflow of the proposed registration method.

2.1. Voxel-based planes generation

In this first step, the source and target point clouds are voxelized with a 3D grid structure and the voxels with planar properties are selected firstly; Then, voxel-based planes will be generated. This step includes two major procedures: (1) planar voxels selection and (2) voxel-based planes calculation.

2.1.1. Planar voxels selection

Step1: We adopt the octree-based voxelization to decompose the entire point cloud with 3D cubic grids. The major advantages of using the octree structure were stated by Boerner, Hoegner, and Stilla (2017) and Xu et al. (2017a). In theory, the smaller is the size of the voxel, the larger is the number of details of the scan that can be preserved after the voxelization, as well as the more is the number of extracted voxel-based planes. In our method, voxels of different point clouds from the same dataset are equal in size, but point clouds from different datasets can be voxelized with different voxel sizes. In our study, the size of the voxel is determined according to applications and manually set, which is one of the two important parameters in the process of the proposed registration method.

Step2: The points in voxel $V_i$ are counted as $V_i$, if

(1) $v_i<{\delta }_n,$(1)

then the voxel will be removed. ${\delta }_n$ is a threshold, which represents the minimum point number in a voxel, and it is another one of the two important parameters. The aim of setting threshold ${\delta }_n$ is to improve the calculation accuracy of the voxel-based planes parameters, because the more points are used to fit a plane, the more robust its parameters are to the noises.

Step3: The voxel property is calculated. Only the voxels with planar property are preserved to reduce the number of the extracted voxel-based planes and improve the accuracy of estimated planar parameters, both of which will highly influent the speed and quality of the registration. The point set inside the voxel $V_i$ is regarded as $\left\{p_1,p_2,\cdots ,p_n\right\}$, then its covariance matrix ${\mathit{M}}_{\mathit{i}}$ can be constructed by Equation (2)(2) ${\mathit{M}}_{\mathit{i}}=\sum _{j=1}^n\left[p_j-\bar{p}\right]{\left[p_j-\bar{p}\right]}^T,$(2)

(2) ${\mathit{M}}_{\mathit{i}}=\sum _{j=1}^n\left[p_j-\bar{p}\right]{\left[p_j-\bar{p}\right]}^T,$(2)

where $p_j$ denotes the point in the voxel; $\bar{p}$ is the centroid of these points; and $n$ is the points number of the voxel $V_i$.

The planar property $c_i$ of voxel $V_i$ is calculated by the eigenvalues from the eigenvalue decomposition (EVD) of the covariance matrix ${\mathit{M}}_{\mathit{i}}$ by Equation (3)(3) $c_i=\frac{{\lambda }_1^i}{{\lambda }_1^i+{\lambda }_2^i+{\lambda }_3^i},$(3)

(3) $c_i=\frac{{\lambda }_1^i}{{\lambda }_1^i+{\lambda }_2^i+{\lambda }_3^i},$(3)

where $\left({\lambda }_1^i,{\lambda }_2^i,{\lambda }_3^i\right)$ in the decreasing order of magnitude represents the eigenvalues from EVD of ${\mathit{M}}_{\mathit{i}}$. If

(4) $c_i<{\delta }_c,$(4)

then the voxel $V_i$ is regarded as a planar voxel and preserved. ${\delta }_c$is a threshold that determines a voxel is a planar one or not, which was set as 0.03 in this paper.

2.1.2. Calculation of voxel-based planes

The Hessian normal form of plane equation is adopted in this work, i.e. a plane is represented as

(5) $\mathit{n}\cdot \mathit{p}-d=0,$(5)

where $\mathit{n}=\left(\begin{array}{ccc}{n}_x& n_y& n_z\end{array}\right)$ is the unit normal vector of the plane; $\mathit{p}={\left(\begin{array}{ccc}x& y& z\end{array}\right)}^T$ is an arbitrary point on the plane; and $d$ is a constant term denoting the perpendicular distance from the origin to the plane.

Regard the plane of the voxel $V_i$ as $P_i$, which four parameters, the unit normal vector $n_i$ and its constant term $d_i$, can be calculated by the random sampling consistency (RANSAC) algorithm using the point set $\left\{p_1,p_2,\cdots ,p_n\right\}$ inside the voxel $V_i$.

Step1: Select three points from the point set $\left\{p_1,p_2,\mathrm{L},p_n\right\}$ randomly, which are used to calculate a plane $P_j$, e.g. ${\mathit{n}}_j\cdot \mathit{p}-d_j=0$.

Step2: Calculate the distance $\mathrm{\Delta }{d}_k\left(k=1,2,\mathrm{L},n\right)$ from the point $p_k$ in $\left\{p_1,p_2,\mathrm{L},p_n\right\}$ to the plane $P_j$. If

(6) $\mathrm{\Delta }{d}_k<{\delta }_d,$(6)

then the point $p_k$ is regarded as inlier; otherwise, as outlier. ${\delta }_d$is a threshold to determine the inliers, and in this paper, it was set as 0.5 m.

Step3: The plane $P_j$ is refined by the inliers obtained in Step 2, and then set the refined plane $P_i$ as the plane of the voxel $V_i$.

Step4: Repeat the above three Steps. Only when the number of the inliers is larger than a threshold, which was set as 0.8 times of the points in the voxel, or the iterations reach the maximum, which was set as 50 in this paper, the iteration is stopped and the final plane $P_i$ is outputted.

In theory, there are two normal vectors that go in opposite directions for any given plane. The choice of normal vector orientation in source and target coordinate systems may result in incorrect solution of the transformation matrix or inconsistent solvable equation. Therefore, the direction of the normal vector needs to be defined in a consistent way (Khoshelham 2016). Here, we just adjust the normal vector to make $d>0$ to eliminate the ambiguity of the direction of the normal vector.

2.2. Conjugate 2-plane base pairs construction

Traditional plane-based registration methods usually select a triple of planes from each scan to form a 3D corner (Xiao et al. 2013; Xu et al. 2017b), which defines a coordinate system that determines six DOFs. By finding the corresponding planes from conjugate triples of planes, transformation parameters can be achieved (Theiler and Schindler 2012). On one hand, there are more geometric constraints between three planes, making it difficult and complex to recognize the conjugate triples of planes from pairwise point clouds; On the other hand, there are six kinds of corresponding plane pairs from two triples of planes, which burdens the calculation of transformation matrix. If one plane was removed from the triple of planes (forming 2-plane base), the geometric constraints between planes and the complexity of plane correspondences constructing (there are only two kinds of corresponding planes between a pair of conjugate 2-plane bases) could be reduced significantly. With this idea, we proposed the congruent voxel-based 2-plane bases for estimating plane correspondences, which are further used to estimate transformation matrix. Compared with the classic triple of planes set base, our method uses neither the simple constraint of the direction of normal vectors nor the positions of intersecting points of planes. Instead, we utilize the invariant congruence formed by the dihedral angle between two planes to find the plane correspondences, which is more efficient and has good correctness.

2.2.1. Generation of 2-plane bases

The plane set $\left\{P_1,P_2,\mathrm{L},P_n\right\}$ extracted from a point cloud is used to generate 2-plane bases via the dihedral angle $\theta$ of two planes. If

(7) ${\delta }_{\theta 1}<{\theta }_{ij}<{\delta }_{\theta 2},$(7)

then the plane $P_i$ and $P_j$ are regarded as a 2-plane base. ${\theta }_{ij}$ represents the dihedral angle between plane $P_i$ and $P_j$; ${\delta }_{\theta 1}$ and ${\delta }_{\theta 2}$ are angle threshold, which were set as 10 degrees and 80 degrees in this paper. Setting ${\delta }_{\theta 1}$ is to avoid to select two approximate parallel planes as a 2-plane base, which cannot be used to calculate the rotation matrix between pairwise point clouds; Setting ${\delta }_{\theta 2}$ is to control the number of generated 2-plane bases, which determines registration efficiency.

As shown in Figure 2(a), each 2-plane base consists of two planes and their dihedral angle $\theta$. In Algorithm 1, we provide a detailed procedure of the generation of the available 2-plane bases in a point cloud.

Table

Algorithm 1: 2-plane bases generation
Input: Voxel-based planes $\left\{P_1,P_2,\cdots ,P_n\right\}$
Output: 2-plane bases $B=\left\{\left(P_i,P_j,{\theta }_{ij}\right),\left(P_h,P_k,{\theta }_{hk}\right),\cdots ,\left(P_m,P_l,{\theta }_{ml}\right)\right\}$
1:	for i = 1 to n-1 do
2:	for j = i+1 to n do
3:	Calculate the dihedral angle ${\theta }_{ij}$ between plane $P_i$ and $P_j$.
4:	if ${\delta }_{\theta 1}<{\theta }_{ij}<{\delta }_{\theta 2}$ do
5:	$\mathit{B}.push\mathrm{\_}back\left(\left(P_i,P_j,{\theta }_{ij}\right)\right)$
6:	end if
7:	end for j
8:	end for i

Figure 2. Definition of 2-plane base and conjugate 2-plane base pair. (a) Definition of 2-plane base; (b) Definition of conjugate 2-plane base pair.

2.2.2. Recognition of initial plane correspondences by conjugate 2-plane base pairs

To identify the initial corresponding voxel-based planes from source and target point clouds, the procedure is divided into three successive steps: (1) selecting conjugate 2-plane base pairs, (2) calculating rotation matrix, and (3) identifying initial plane correspondences.

• Selecting conjugate 2-plane base pairs

As shown in Figure 2(b), for each 2-plane base $\left({\hspace{0.17em}}^s{P}_i^{\hspace{0.17em}},{\hspace{0.17em}}^s{P}_j^{\hspace{0.17em}},{\hspace{0.17em}}^s{\theta }_{ij}^{\hspace{0.17em}}\right)$ in source point cloud, its nearest neighbour $\left({\hspace{0.17em}}^t{P}_h^{\hspace{0.17em}},{\hspace{0.17em}}^t{P}_k^{\hspace{0.17em}},{\hspace{0.17em}}^t{\theta }_{hk}^{\hspace{0.17em}}\right)$ from target point cloud is found by comparing ${\hspace{0.17em}}^s{\theta }_{ij}^{\hspace{0.17em}}$ and${\hspace{0.17em}}^t{\theta }_{hk}^{\hspace{0.17em}}$. Assuming that the $\left({\hspace{0.17em}}^s{P}_i^{\hspace{0.17em}},{\hspace{0.17em}}^s{P}_j^{\hspace{0.17em}},{\hspace{0.17em}}^s{\theta }_{ij}^{\hspace{0.17em}}\right)$ and $\left({\hspace{0.17em}}^t{P}_h^{\hspace{0.17em}},{\hspace{0.17em}}^t{P}_k^{\hspace{0.17em}},{\hspace{0.17em}}^t{\theta }_{hk}^{\hspace{0.17em}}\right)$ are a pair of conjugate 2-plane bases, the plane correspondences can be $\left({\hspace{0.17em}}^s{P}_i^{\hspace{0.17em}},{\hspace{0.17em}}^s{P}_j^{\hspace{0.17em}}\right)$ and $\left({\hspace{0.17em}}^t{P}_h^{\hspace{0.17em}},{\hspace{0.17em}}^t{P}_k^{\hspace{0.17em}}\right)$, but they can also be $\left({\hspace{0.17em}}^s{P}_i^{\hspace{0.17em}},{\hspace{0.17em}}^s{P}_j^{\hspace{0.17em}}\right)$ and $\left({\hspace{0.17em}}^{\mathrm{t}}{P}_k^{\hspace{0.17em}},{\hspace{0.17em}}^{\mathrm{t}}{P}_h^{\hspace{0.17em}}\right)$. Therefore, we can calculate two rotation matrixes using a pair of conjugate 2-plane bases.

• Calculating rotation matrix

Taking $\left({\hspace{0.17em}}^s{P}_i^{\hspace{0.17em}},{\hspace{0.17em}}^s{P}_j^{\hspace{0.17em}}\right)$ and $\left({\hspace{0.17em}}^t{P}_h^{\hspace{0.17em}},{\hspace{0.17em}}^t{P}_k^{\hspace{0.17em}}\right)$ as example, we can formulate

(8) $\mathit{R}\left(\begin{array}{c}\begin{array}{cc}{\hspace{0.17em}}^s{n}_x^i& {\hspace{0.17em}}^s{n}_x^j\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^s{n}_y^i& {\hspace{0.17em}}^s{n}_y^j\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^s{n}_z^i& {\hspace{0.17em}}^s{n}_z^j\end{array}\end{array}\right)=\left(\begin{array}{c}\begin{array}{cc}{\hspace{0.17em}}^t{n}_x^h& {\hspace{0.17em}}^t{n}_x^k\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^t{n}_y^h& {\hspace{0.17em}}^t{n}_y^k\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^t{n}_z^h& {\hspace{0.17em}}^t{n}_z^k\end{array}\end{array}\right),$(8)

where $R$ represents the 3D rotation matrix between the source and target point clouds; and $\left(\begin{array}{ccc}{n}_x^l& n_y^l& n_z^l\end{array}\right)$ represents the unit normal vector of the plane $P_l$.

Performing the singular value decomposition (SVD) for the Equation (8)(8) $\mathit{R}\left(\begin{array}{c}\begin{array}{cc}{\hspace{0.17em}}^s{n}_x^i& {\hspace{0.17em}}^s{n}_x^j\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^s{n}_y^i& {\hspace{0.17em}}^s{n}_y^j\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^s{n}_z^i& {\hspace{0.17em}}^s{n}_z^j\end{array}\end{array}\right)=\left(\begin{array}{c}\begin{array}{cc}{\hspace{0.17em}}^t{n}_x^h& {\hspace{0.17em}}^t{n}_x^k\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^t{n}_y^h& {\hspace{0.17em}}^t{n}_y^k\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^t{n}_z^h& {\hspace{0.17em}}^t{n}_z^k\end{array}\end{array}\right),$(8)

(9) $\left[\begin{array}{ccc}\mathit{U}& \mathit{S}& \mathit{V}\end{array}\right]=\mathrm{s}\mathrm{v}\mathrm{d}\left(\left(\begin{array}{c}\begin{array}{cc}{\hspace{0.17em}}^s{n}_x^i& {\hspace{0.17em}}^s{n}_x^j\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^s{n}_y^i& {\hspace{0.17em}}^s{n}_y^j\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^s{n}_z^i& {\hspace{0.17em}}^s{n}_z^j\end{array}\end{array}\right){\left(\begin{array}{c}\begin{array}{cc}{\hspace{0.17em}}^t{n}_x^h& {\hspace{0.17em}}^t{n}_x^k\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^t{n}_y^h& {\hspace{0.17em}}^t{n}_y^k\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^t{n}_z^h& {\hspace{0.17em}}^t{n}_z^k\end{array}\end{array}\right)}^T\right),$(9)

then the 3D rotation matrix $\mathit{R}$ can be represented as

(10) $\mathit{R}=\mathit{V}{\mathit{U}}^T.$(10)

Using a rotation matrix $\mathit{R}$ that calculated by a pair of conjugate 2-plane bases, a set of initial plane correspondences can be identified.

• Identifying initial plane correspondences

The corresponding voxel-base planes from the pairwise point clouds are identified using the unit normal vector of the plane. Assuming the original unit normal vectors of the planes in source and target point clouds as ${\hspace{0.17em}}^{s}n=\left({\hspace{0.17em}}^s{\mathit{n}}_1^{\hspace{0.17em}},{\hspace{0.17em}}^s{\mathit{n}}_2^{\hspace{0.17em}},\mathrm{L},{\hspace{0.17em}}^s{\mathit{n}}_{m1}^{\hspace{0.17em}}\right)$ and ${\hspace{0.17em}}^t\mathit{n}=\left({\hspace{0.17em}}^t{\mathit{n}}_1^{\hspace{0.17em}},{\hspace{0.17em}}^t{\mathit{n}}_2^{\hspace{0.17em}},\mathrm{L},{\hspace{0.17em}}^t{\mathit{n}}_{m2}^{\hspace{0.17em}}\right)$ respectively, where $n_i=\left(\begin{array}{ccc}{n}_x& n_y& n_z\end{array}\right)$ represents the unit normal vector of the i-th plane, and m1 and m2 are the number of the planes in source and target point clouds, respectively. Then, the initial plane correspondences are identified using ${\hspace{0.17em}}^{s}n$, ${\hspace{0.17em}}^{t}n$ and $R$.

Step1: The planes in source point cloud are rotated using the rotation matrix $R$, and the unit normal vector of the rotated planes are regarded as ${\hspace{0.17em}}_{\mathrm{R}}^{s}n=\left({\hspace{0.17em}}_{\mathrm{R}}^s{n}_1^{\hspace{0.17em}},{\hspace{0.17em}}_{\mathrm{R}}^s{n}_2^{\hspace{0.17em}},\mathrm{L},{\hspace{0.17em}}_{\mathrm{R}}^s{n}_{m1}^{\hspace{0.17em}}\right)$.

Step2: For each plane normal vector ${\hspace{0.17em}}_{\mathrm{R}}^s{n}_h^{\hspace{0.17em}}$ in ${\hspace{0.17em}}_{\mathrm{R}}^{s}n$, if there exists an element ${\hspace{0.17em}}^t{n}_k^{\hspace{0.17em}}$ in ${\hspace{0.17em}}^{t}n$ that satisfies the constraint as in Equation (11)(11) $\left\{\begin{array}{c}k=argmin∥{\hspace{0.17em}}_R^s{\mathit{n}}_h^{\hspace{0.17em}}-{\hspace{0.17em}}^t{\mathit{n}}_i^{\hspace{0.17em}}∥\left(i=1,2,\mathrm{L},m2\right)\\ h=argmin∥{\hspace{0.17em}}_R^s{\mathit{n}}_i^{\hspace{0.17em}}-{\hspace{0.17em}}^t{\mathit{n}}_k^{\hspace{0.17em}}∥\left(i=1,2,\mathrm{L},m1\right)\end{array}\right.,$(11) , then the corresponding planes ${\hspace{0.17em}}^s{P}_h^{\hspace{0.17em}}$ and ${\hspace{0.17em}}^t{P}_k^{\hspace{0.17em}}$ are considered as a pair of corresponding planes. In other words, only when ${\hspace{0.17em}}^t{n}_k^{\hspace{0.17em}}$ was the nearest neighbour to ${\hspace{0.17em}}_{\mathrm{R}}^s{n}_h^{\hspace{0.17em}}$ in ${\hspace{0.17em}}^{t}n$, and ${\hspace{0.17em}}_{\mathrm{R}}^s{n}_h^{\hspace{0.17em}}$ to ${\hspace{0.17em}}^t{n}_k^{\hspace{0.17em}}$ in ${\hspace{0.17em}}_{\mathrm{R}}^{s}n$, then ${\hspace{0.17em}}^s{P}_h^{\hspace{0.17em}}$ and ${\hspace{0.17em}}^t{P}_k^{\hspace{0.17em}}$ were identified as an initial plane correspondence. All the initial plane correspondences are regarded as $C=\left\{C_1,C_2,\mathrm{L},C_m\right\}$, where $C_i=\left\{{\hspace{0.17em}}^s{P}_h^{\hspace{0.17em}},{\hspace{0.17em}}^t{P}_k^{\hspace{0.17em}}\right\}$ is the i-th initial plane correspondence, and m represents the number of the initial plane correspondences.

(11) $\left\{\begin{array}{c}k=argmin∥{\hspace{0.17em}}_R^s{\mathit{n}}_h^{\hspace{0.17em}}-{\hspace{0.17em}}^t{\mathit{n}}_i^{\hspace{0.17em}}∥\left(i=1,2,\mathrm{L},m2\right)\\ h=argmin∥{\hspace{0.17em}}_R^s{\mathit{n}}_i^{\hspace{0.17em}}-{\hspace{0.17em}}^t{\mathit{n}}_k^{\hspace{0.17em}}∥\left(i=1,2,\mathrm{L},m1\right)\end{array}\right.,$(11)

where the norm is the 3D Euclidean distance.

By more than three plane correspondences, the transformation matrix of source and target point clouds can be achieved using the parameters of the corresponding planes, e.g. the unit normal vector $\mathit{n}=\left(\begin{array}{ccc}{n}_x& n_y& n_z\end{array}\right)$ and constant term d.

2.3. Registration using corresponding planes

The registration consists of two essential steps: (1) estimation of initial transformation parameters and (2) optimization of final transformation parameters. The initial transformation parameters are estimated by the initial plane correspondences identified in Section 2.2. The optimized registration parameters are estimated using optimal corresponding planes recognized by the initial transformation matrix.

2.3.1. Estimation of initial transformation parameters

With the initial plane correspondences $C=\left\{C_1,C_2,\mathrm{L},C_m\right\}$, the following Equation can be formulated

(12) $\left\{\begin{array}{c}{{\mathit{R}}_{ini}}^{\mathrm{s}}\mathit{n}={\hspace{0.17em}}^t\mathit{n}\\ {\hspace{0.17em}}^t{\mathit{n}}^T{\mathit{t}}_{ini}={\left[{\hspace{0.17em}}^s\mathit{d}-{\hspace{0.17em}}^t\mathit{d}\right]}^T\end{array}\right.,$(12)

where ${ }^{s}n=\left({ }^s{\mathit{n}}_1^{\hspace{0.17em}}{,}^s{\mathit{n}}_2^{\hspace{0.17em}},\mathrm{L}{,}^s{\mathit{n}}_m^{\hspace{0.17em}}\right)$ and ${ }^{t}n=\left({ }^t{\mathit{n}}_1^{\hspace{0.17em}}{,}^t{\mathit{n}}_2^{\hspace{0.17em}},\mathrm{L}{,}^t{\mathit{n}}_m^{\hspace{0.17em}}\right)$ are the unit normal vectors of the plane from the source and target point clouds, $n_i=\left(\begin{array}{ccc}{n}_x& n_y& n_z\end{array}\right)$ represents the unit normal vector of the i-th plane; ${ }^{\mathrm{s}}\mathit{d}=\left({ }^{\mathrm{s}}{d}_1^{\hspace{0.17em}}{,}^{\mathrm{s}}{d}_2^{\hspace{0.17em}},\mathrm{L}{,}^{\mathrm{s}}{d}_m^{\hspace{0.17em}}\right)$ and ${ }^t\mathit{d}=\left({ }^t{d}_1^{\hspace{0.17em}}{,}^t{d}_2^{\hspace{0.17em}},\mathrm{L}{,}^t{d}_m^{\hspace{0.17em}}\right)$ are the constant term; ${\mathit{R}}_{ini}$ and ${\mathit{t}}_{ini}$ are the rotation and translation component of the initial transformation matrix.

The rotation component ${\mathit{R}}_{ini}$ can be obtained by the SVD method using Equation (8)(8) $\mathit{R}\left(\begin{array}{c}\begin{array}{cc}{\hspace{0.17em}}^s{n}_x^i& {\hspace{0.17em}}^s{n}_x^j\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^s{n}_y^i& {\hspace{0.17em}}^s{n}_y^j\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^s{n}_z^i& {\hspace{0.17em}}^s{n}_z^j\end{array}\end{array}\right)=\left(\begin{array}{c}\begin{array}{cc}{\hspace{0.17em}}^t{n}_x^h& {\hspace{0.17em}}^t{n}_x^k\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^t{n}_y^h& {\hspace{0.17em}}^t{n}_y^k\end{array}\\ \begin{array}{cc}{\hspace{0.17em}}^t{n}_z^h& {\hspace{0.17em}}^t{n}_z^k\end{array}\end{array}\right),$(8) -(10(10) $\mathit{R}=\mathit{V}{\mathit{U}}^T.$(10) ), and the translation component ${\mathit{t}}_{ini}$ can be solved using the least square (LS) algorithm, as shown in Equation (13)(13) ${\mathit{t}}_{ini}={\left({\hspace{0.17em}}^t\mathit{n}{\hspace{0.17em}}^t{\mathit{n}}^T\right)}^{-1}{\hspace{0.17em}}^t\mathit{n}\left({\left[{\hspace{0.17em}}^s\mathit{d}-{\hspace{0.17em}}^t\mathit{d}\right]}^T\right).$(13)

(13) ${\mathit{t}}_{ini}={\left({\hspace{0.17em}}^t\mathit{n}{\hspace{0.17em}}^t{\mathit{n}}^T\right)}^{-1}{\hspace{0.17em}}^t\mathit{n}\left({\left[{\hspace{0.17em}}^s\mathit{d}-{\hspace{0.17em}}^t\mathit{d}\right]}^T\right).$(13)

Combing $R_{ini}$ and $t_{ini}$, the initial transform matrix $T_{ini}$ can be obtained

(14) ${\mathit{T}}_{ini}=\left[\begin{array}{cc}\begin{array}{c}{\mathit{R}}_{ini}\\ 0\end{array}& \begin{array}{c}{\mathit{t}}_{ini}\\ 1\end{array}\end{array}\right],$(14)

then its corresponding LCP can be calculated and optimal plane correspondences can be obtained based on $T_{ini}$.

Step1: Assuming the plane parameters as $\pi =[\begin{array}{cccc}{n}_x& n_y& n_z& d\end{array}]$, and the planes in the source point cloud are transformed by ${\mathit{T}}_{ini}$

(15) ${ }_T^{\text{s}}{\pi }^{ }{=}_{ }^{\text{s}}\pi T_{ini}.$(15)

where ${\hspace{0.17em}}_{\hspace{0.17em}}^s\pi$ and ${\hspace{0.17em}}_T^s\pi$ represent the original planes and the transformed ones from the source point cloud, respectively.

Step2: Based on the initial transformation matrix $T_{ini}$, the optimal plane correspondences are identified. First, the nearest planes are calculated from ${\hspace{0.17em}}_{\mathrm{T}}^{\mathrm{s}}\pi$ and ${\hspace{0.17em}}^{\mathrm{t}}\pi$ (${\hspace{0.17em}}^{\mathrm{t}}\pi$ representing the planes of target point cloud) based on their unit normal vector by the Equation (11)(11) $\left\{\begin{array}{c}k=argmin∥{\hspace{0.17em}}_R^s{\mathit{n}}_h^{\hspace{0.17em}}-{\hspace{0.17em}}^t{\mathit{n}}_i^{\hspace{0.17em}}∥\left(i=1,2,\mathrm{L},m2\right)\\ h=argmin∥{\hspace{0.17em}}_R^s{\mathit{n}}_i^{\hspace{0.17em}}-{\hspace{0.17em}}^t{\mathit{n}}_k^{\hspace{0.17em}}∥\left(i=1,2,\mathrm{L},m1\right)\end{array}\right.,$(11) ; Then, the constant term of the plane parameters is used to confirm the nearest planes are corresponding planes or not. If the nearest planes satisfy the constraint in Equation (16)(16) $∥{\hspace{0.17em}}_T^{\mathrm{s}}{d}_h^{\hspace{0.17em}}-{\hspace{0.17em}}^{\mathrm{t}}{d}_k^{\hspace{0.17em}}∥<{\delta }_m,$(16) , then they are regarded as a pair of corresponding planes, otherwise not. The values of LCP corresponding to the initial transformation matrix is just equal to the number of the plane correspondences, and all the corresponding planes form the set of optimal plane correspondences.

(16) $∥{\hspace{0.17em}}_T^{\mathrm{s}}{d}_h^{\hspace{0.17em}}-{\hspace{0.17em}}^{\mathrm{t}}{d}_k^{\hspace{0.17em}}∥<{\delta }_m,$(16)

where ${\delta }_m$ is a threshold, which was set as 1 m in this paper.

2.3.2. Recognition of initial plane correspondences by conjugate 2-plane base pairs

Although we have computed the initial transformation parameters $R_{ini}$ and $t_{ini}$, it still needs to estimate the optimal transformation parameters $R_{opt}$ and $t_{opt}$. Here, we adopt the method of applying the optimal plane correspondences to get the closed-form solution of transformation matrix estimating (Khoshelham 2016; Forstner and Khoshelham 2017). It is worth mentioning that for our closed-form solution, we do not apply the normalization using mean-centring and scaling, because translation normalization (i.e. mean-centring) has been realized in the process of initial translation parameter estimation, and scaling is not considered in our study, because the LiDAR technology obtains the real distance between points.

Based on the set of optimal plane correspondences, the Equation (12)(12) $\left\{\begin{array}{c}{{\mathit{R}}_{ini}}^{\mathrm{s}}\mathit{n}={\hspace{0.17em}}^t\mathit{n}\\ {\hspace{0.17em}}^t{\mathit{n}}^T{\mathit{t}}_{ini}={\left[{\hspace{0.17em}}^s\mathit{d}-{\hspace{0.17em}}^t\mathit{d}\right]}^T\end{array}\right.,$(12) can be updated, and the optimum transformation parameters $R_{opt}$ and $t_{opt}$ can be calculated by the SVD and LS algorithm, respectively.

By the before-mentioned steps of the proposed coarse registration scheme, we know that with a pair of conjugate 2-plane bases, an initial transformation matrix as well as its LCP can be calculated. For each 2-plane base in the source point cloud, its corresponding 2-plane base can be found from the target point cloud using the traversing strategy; Then a series of initial transformation matrixes with LCPs are obtained; Finally, the transformation matrix with the highest LCP is optimized using its generated optimal plane correspondences, and the optimized transformation matrix is outputted as the final coarse registration result. To illustrate our method more clearly, in Algorithm 2, we provide a detailed procedure for the proposed coarse registration algorithm.

Table

Algorithm 2: Pairwise coarse registration by traversing voxel-based 2-plane bases
Input: Source and target point clouds, $\left\{{\hspace{0.17em}}^{\mathrm{s}}{p}_1^{\hspace{0.17em}},{\hspace{0.17em}}^{\mathrm{s}}{p}_2^{\hspace{0.17em}},\cdots ,{\hspace{0.17em}}^{\mathrm{s}}{p}_{n1}^{\hspace{0.17em}}\right\}$ and $\left\{{\hspace{0.17em}}^{\mathrm{t}}{p}_1^{\hspace{0.17em}},{\hspace{0.17em}}^{\mathrm{t}}{p}_2^{\hspace{0.17em}},\cdots ,{\hspace{0.17em}}^{\mathrm{t}}{p}_{n2}^{\hspace{0.17em}}\right\}$.
Output: The coarse registration matrix $T$.
1:	Extract the voxel-based planes from source and target point clouds, regarded as $\left\{{\hspace{0.17em}}^{\mathrm{s}}{P}_1^{\hspace{0.17em}},{\hspace{0.17em}}^{\mathrm{s}}{P}_2^{\hspace{0.17em}},\cdots ,{\hspace{0.17em}}^{\mathrm{s}}{P}_{m1}^{\hspace{0.17em}}\right\}$ and $\left\{{\hspace{0.17em}}^{\mathrm{t}}{P}_1^{\hspace{0.17em}},{\hspace{0.17em}}^{\mathrm{t}}{P}_2^{\hspace{0.17em}},\cdots ,{\hspace{0.17em}}^{\mathrm{t}}{P}_{m2}^{\hspace{0.17em}}\right\}$.
2:	Generation the 2-plane bases from the voxel-based planes in source and target point clouds, regarded as $\left\{{\hspace{0.17em}}^{\mathrm{s}}{A}_1^{\hspace{0.17em}},{\hspace{0.17em}}^{\mathrm{s}}{A}_2^{\hspace{0.17em}},\cdots ,{\hspace{0.17em}}^{\mathrm{s}}{A}_{k1}^{\hspace{0.17em}}\right\}$ and $\left\{{\hspace{0.17em}}^{\mathrm{t}}{A}_1^{\hspace{0.17em}},{\hspace{0.17em}}^{\mathrm{t}}{A}_2^{\hspace{0.17em}},\cdots ,{\hspace{0.17em}}^{\mathrm{t}}{A}_{k2}^{\hspace{0.17em}}\right\}$.
3:	for i = 1 to k1 do
4:	for j = 1 to k2 do
5:	Find a corresponding 2-plane base for ${\hspace{0.17em}}^{\mathrm{s}}{A}_i^{\hspace{0.17em}}$ from $\left\{{\hspace{0.17em}}^{\mathrm{t}}{A}_1^{\hspace{0.17em}},{\hspace{0.17em}}^{\mathrm{t}}{A}_2^{\hspace{0.17em}},\cdots ,{\hspace{0.17em}}^{\mathrm{t}}{A}_{k2}^{\hspace{0.17em}}\right\}$.
6:	Calculate the rotation matrix R by the pair of conjugate 2-plane bases.
7:	Find the initial plane correspondences $C_{ini}$ by R.
8:	Calculate the initial transformation matrix $T_{ini}$ and its $LC{P}_i$ using $C_{ini}$.
9:	Find the optimal plane correspondences $C_{opt}$ by $T_{ini}$.
10:	if $LC{P}_i>LC{P}_0$ do
11:	$LC{P}_0=LC{P}_i$.
12:	end if
13:	Obtain the optimized coarse transformation matrix $T_{\text{user}2opt}$ using $C_{\text{user}2opt}$.
14:	end for j
15:	end for i

3. Experimental results and discussion

In this section, the experiments are performed on five datasets, two of which belong to the indoor dataset, and the other three are the outdoor dataset, to test the result of the proposed method in terms of time performance, registration accuracy, root mean square distance (RMSD) and successful registration rate (SRR), as well as comparison with other state-of-the-art algorithms. All the experiments were carried out on a computer with Windows 10, Intel(R) Xeon(R) E-2176 G CPU at 3.70 GHz and 64.0-GB RAM.

3.1. Experimental datasets and evaluation metric

3.1.1. Experimental datasets

The performance of the proposed registration algorithm is evaluated using five datasets (i.e. RESSO(6i), RESSO(6j), RESSO(7d), RESSO(7e), and Bremen city), which include ground-truths that allows for a proper assessment. The RESSO benchmark is a TLS dataset presented by Chen et al. (2020) to evaluate registering point clouds with a small overlap. The Bremen city benchmark is dataset from the Thermal Mapper project that was acquired by the Jacobs University Bremen (Borrmann, Elseberg, and Nüchter 2013). Table 1 provides a detailed description of the five datasets that used in this paper, and Figure 3 shows the full view of the experimental datasets.

Figure 3. Full view of the five experimental datasets, with shading representing differences in height. (a), (b), (c), (d) and (e) show the RESSO(6i), RESSO(6j), RESSO(7d), RESSO(7e) and Bremen city datasets respectively.

Table 1. Detailed information of the experimental datasets.

Datasets	#Scans	#Points (10×6)	Dimensions (m)
RESSO(6i)	5	1.1	43 × 83 × 19
RESSO(6j)	9	10.7	163 × 243 × 37
RESSO(7d)	3	0.7	338 × 244 × 153
RESSO(7e)	26	4.2	549 × 562 × 84
Bremen city	13	11.4	816 × 841 × 154

3.1.2. Experimental metric

The registration performance is evaluated in terms of registration accuracy, RMSD, SRR, and time performance. The above aspects are commonly adopted by many researchers (Ghorbani et al. 2022; Fu et al. 2021a, 2021b; Xu et al. 2019).

The registration accuracy includes rotation error $e_r$ and translation error $e_t$, which can be calculated by the comparison with the ground-truth. Given the target point cloud $Q$, the transformation matrix $T$ from the source point cloud $P$ to $Q$ can be computed via the registration method. The residual transformation $\mathrm{\Delta }T$from $P$ to $Q$ is defined as

(17) $\mathrm{\Delta }\mathit{T}=\mathit{T}{\left({\mathit{T}}_G\right)}^{-1}=\left[\begin{array}{cc}\mathrm{\Delta }\mathit{R}& \mathrm{\Delta }\mathit{t}\\ 0& 1\end{array}\right],$(17)

where $T_G$ is the ground-truth.

The rotation error $e_r$ and translation error $e_t$ can be calculated based on their corresponding rotational component $\mathrm{\Delta }R$ and translational component$\mathrm{\Delta }t$, as expressed in Equation (18)(18) $\left\{\begin{array}{c}{e}_r=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{trace\left(\mathrm{\Delta }\mathit{R}\right)-1}{2}\right)\\ e_t=∥\mathrm{\Delta }\mathit{t}∥\end{array}\right.,$(18)

(18) $\left\{\begin{array}{c}{e}_r=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{trace\left(\mathrm{\Delta }\mathit{R}\right)-1}{2}\right)\\ e_t=∥\mathrm{\Delta }\mathit{t}∥\end{array}\right.,$(18)

where $trace\left(\mathrm{\Delta }R\right)$ denotes the trace of $\mathrm{\Delta }R$.

The RMSD is a statistical criterion, which is calculated between the source point cloud after applying the registration algorithm (taken as $P^{\text{user}}2T$) and the same point cloud in its true position (taken as $P^{\text{user}}2G$). The RMSD is computed as follows

(19) $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{D}=\frac{\sum _{i=0}^n∥p_i^{\mathrm{T}}-p_i^{\mathrm{G}}∥}{n},$(19)

where $∥\cdot ∥$ is Euclidean distance and n is the number of source point cloud.

The SRR is used as the percentage of completeness, which indicates how many of the point cloud registration was successful and whether it is possible to enter the fine registration stage. This criterion is defined as expressed in Equation (20)(20) $\mathrm{S}\mathrm{R}\mathrm{R}=\frac{N_S}{N}\times 100,$(20)

(20) $\mathrm{S}\mathrm{R}\mathrm{R}=\frac{N_S}{N}\times 100,$(20)

where $N_S$ and $N$ are the number of successful registration and the number of two adjacent point clouds for registration, respectively. (In this paper, a registration result is regarded as a successful one, only when its RMSD is less than 1.0 m).

3.2. Registration result

To clearly demonstrate the registration procedure of the proposed method, in Figure 4, we show the results of extracting voxel-based planes, generating 2-plane bases, calculating the rotation matrix by a pair of conjugate 2-plane bases, and computing the transformation matrix using the conjugate voxel-based planes that are recognized by the obtained rotation matrix.

Figure 4. The registration results of the proposed method in each step. (a) Original input point clouds; (b) Voxel-based planes; (c) 2-plane bases; (d) Aligned point clouds by the rotation matrix thar are calculated by the conjugate 2-plane base pairs; (e) Aligned point clouds by the transformation matrix that are calculated by the conjugate voxel-based planes.

In Figure 4(a), it is the input point clouds. Figure 4(b) shows the voxel-based planes extracting results. Then, a pair of conjugate 2-plane bases is recognized from the pairwise point clouds, as shown in Figure 4(c). In Figure 4(d), it shows the aligned point clouds using the rotation matrix that calculated by a pair of conjugate 2-plane bases, where we can see that it consisted between the directions of the aligned pairwise point clouds, but there are obvious displacements along the x-axis, y-axis and z-axis. The reason is that a pair of conjugate 2-plane bases with five necessary observations (conditions) can only calculate a rotation matrix with three DOFs for a pair of point clouds, it cannot compute a transformation matrix with six DOFs for a pair of point clouds. In other words, a pair of conjugate 2-plane bases can only compute the rotation matrix, but not the translation vector between a pair of point clouds. Therefore, more conjugate voxel-based planes (at least three) are recognized according to the obtained rotation matrix, which are used to perfectly align the pairwise point clouds, as shown in Figure 4(e).

Figure 5 shows a few pairs of point clouds registered by our method on the five different experimental datasets. As we can see from Figure 5, the pairwise point clouds are automatically aligned. There is no significant inconsistency between the two registered point clouds; therefore, it can be qualitatively concluded that the proposed method can feasibly register pairwise point clouds.

Figure 5. Pairs of point clouds registered by our method. (First Two Columns) Input scan pairs. (Right Column) Registration results.

3.3. Evaluation and comparison

To quantitatively evaluate the proposed method for point clouds registration, the registration results are compared with ground-truths in terms of rotation and translation errors, RMSD, SRR, and time performance. To further analyse the performance of our method, two well-known initial alignment algorithms are introduced as baseline for performance comparison, namely, the point-based method—-FPFH (Rusu, Blodow, and Beetz 2009) and the line-based method—-2DLINE (Tao et al. 2020b).

3.3.1. Rotation and translation error evaluation

Table 2 summarizes the rotation error and translation error (i.e., minimum, maximum, average, and the Root Mean Square Error (RMSE)) on the five testing datasets. Note that, only the successful registrations are utilized to calculate the average and RMSE of rotation and translation errors (It is the same for RMSD in the following calculation). These experimental results demonstrate that the proposed method performs well in registering pairwise point clouds, with average rotation error less than 0.4° and translation error less than 0.4 m, which can satisfy the requirements of object extraction and 3D reconstruction.

Table 2. Quantitative evaluation of the proposed method by rotation and translation error.

Datasets	Rotation error (°)				Translation error (m)
	Min	Max	Ave	RMSE	Min	Max	Ave	RMSE
RESSO(6i)	0.16	0.20	0.18	0.02	0.103	1.503	0.206	0.135
RESSO(6j)	0.01	0.28	0.13	0.12	0.047	0.276	0.144	0.079
RESSO(7d)	0.20	0.24	0.22	0.02	0.041	0.052	0.046	0.005
RESSO(7e)	0.01	0.42	0.12	0.11	0.015	0.440	0.157	0.108
Bremen city	0.10	1.02	0.32	0.25	0.063	0.657	0.345	0.183

3.3.2. SRR, RMSD evaluation and time performance

Table 3 reports the SRR, RMSD and time performance of the proposed method, where we can find that:

1. Our method can perfectly register the pairs of point clouds from RESSO(6j), RESSO(7d), RESSO(7e) and Bremen city datasets, but obtain a failing one from RESSO(6i) dataset. The total SRR of the proposed method is about 98%.

2. The average RMSD is all less than 0.42 m, which demonstrates the good performance of our method.

3. It takes only a few minutes to register millions of points, which shows that the proposed method acquires high efficiency in registering large-scale and large amounts of point clouds.

Table 3. Quantitative evaluation of the proposed method by SRR, RMSD and time.

Datasets	SRR (%)	RMSD (m)				Time (s)
		Min	Max	Ave	RMSE
RESSO(6i)	75	0.105	1.501	0.206	0.136	208.9
RESSO(6j)	100	0.503	0.228	0.124	0.065	279.4
RESSO(7d)	100	0.085	0.098	0.091	0.006	62.3
RESSO(7e)	100	0.037	0.324	0.137	0.072	54.6
Bremen city	100	0.105	0.856	0.415	0.205	210.3

In a word, the proposed method can well align pairwise point clouds obtained from indoor (RESSO(6i) and RESSO(6j) datasets) and outdoor (RESSO(7d), RESSO(7e) and Bremen city datasets) scenes.

3.3.3. Performance comparison and analysis

The point-based method—-FPFH (Rusu, Blodow, and Beetz 2009) and the line-based method—-2DLINE (Tao et al. 2020b) are selected as the benchmark methods to conduct the comparing experiments. Table 4 lists their SRR, rotation and translation errors, RMSD and the runtime. It is found that the proposed method outperformed the state-of-the-art registration methods on the testing datasets. More specifically, a number of observations can be noted based on these comparison results.

1. The proposed method (the voxel-plane-based method) outperformed the line-based method (2DLINE method (Tao et al. 2020b)) and the point-based method (FPFH method (Rusu, Blodow, and Beetz 2009)) in terms of SRR, registration errors and RMSD, which demonstrates that voxel-based planes are better elements of registration than feature lines and feature points.

2. Our method achieved the best efficiency, which is attributable to the following two factors. First, the number of generated voxel-based planes is less than the feature lines generated by 2DLINE method and the feature points used in FPFH method. Second, compared with FPFH method, which spends a lot of time to calculate the distinctive descriptors for the feature points using their neighbours, the proposed method just needs to decompose the entire point cloud with 3D cubic grids by the octree-based voxelization and calculate the voxel-based planes, which is a more time-efficiency procedure.

Table 4. Performance comparison.

Datasets	Method	SRR (%)	Rotation error (°)		Translation error (m)		RMSD (m)		Time (s)
			Ave	RMSE	Ave	RMSE	Ave	RMSE
RESSO(6i)	FPFH	50	45.1	57.36	3.613	3.498	5.096	4.221	342.6
	2DLINE	75	1.66	0.75	2.858	3.646	2.839	3.639	194.6
	Ours	75	0.18	0.02	0.206	0.135	0.206	0.136	208.9
RESSO(6j)	FPFH	17	115.40	70.93	25.105	11.533	28.913	13.138	1691.1
	2DLINE	17	89.99	51.88	14.913	10.276	27.841	13.837	8616.3
	Ours	100	0.13	0.12	0.144	0.079	0.124	0.065	279.4
RESSO(7d)	FPFH	50	5.34	1.65	2.074	0.623	3.003	1.001	305.4
	2DLINE	100	0.78	0.22	0.505	0.423	0.610	0.415	159.7
	Ours	100	0.22	0.02	0.046	0.005	0.091	0.006	62.3
RESSO(7e)	FPFH	32	3.45	1.85	3.737	2.419	2.721	1.445	259.6
	2DLINE	92	0.23	0.61	0.410	1.071	0.379	0.727	348.4
	Ours	100	0.12	0.11	0.157	0.108	0.137	0.072	54.6
Bremen city	FPFH	25	2.29	1.25	3.543	1.443	4.260	1.907	1311.5
	2DLINE	83	2.38	1.10	2.433	1.852	1.243	1.115	10571.1
	Ours	100	0.32	0.25	0.345	0.183	0.415	0.205	210.3

4. Conclusions

In this paper, we proposed a target-less algorithm for fast and coarse registration of pairwise point clouds. The voxel-based planes are separated from the original point cloud using the octree-based voxelization, then the 2-plane bases are constructed. After that, the conjugate 2-plane bases are generated by the dihedral angle from the pairwise point clouds and used to calculate the initial plane correspondences, which are further used to calculate the initial transformation matrix as well as the corresponding LCP. The matrix with the highest LCP is applied to find the optimal plane correspondences, which are further utilized to optimize the former transformation matrix. Based on the geometric feature of points within the voxel structure, the number of selected plane feature is significantly less than the point number, which can improve the efficiency of registration processing compared with traditional point-based methods. And the use of plane feature as basic element for registration can be more robust when handling point clouds with varying point density. Compared with line-based method, which lacks the theoretical support because it cannot determine the corresponding points by the line feature for calculating the translation vector, the proposed registration method using the voxel-based plane as basic element has strict support of mathematical theory. The experimental results using different point clouds reveal that our proposed method is effective and outperforms the baseline methods. For the task of fast orientation between two scans, our proposed method can achieve rotation error less than 0.4 degrees, translation error less 0.4 m, RMSD less than 0.42 m, and SRR about 98%. The registration results satisfy the goal of coarse registration. In the future, we will study on the automatic selecting of the voxel size and minimum points, and add robust estimation to the process of translation vector calculating to further improve its accuracy.

Author contributions

Z. Li. designed the method. Y. Fu. realized the method and wrote the original paper. F. Xiong. got the results of the experiments. H. He. analysed these results. W. Wang. and Y. Deng. checked and revised the paper. All the authors have read and agreed to publish this version of the manuscript.

Acknowledgments

We would like to thank the people of the algorithms and open source code involved in this paper, especially Wuyong Tao for sharing his ideas about line-based registration algorithm and code to us. Meanwhile, we would also like to thank the providers of the datasets used in this paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).