An Effective Approach for Point Clouds Registration Based on the Hard and Soft Assignments

TL;DR

This paper introduces a novel point cloud registration method that combines hard and soft assignment strategies to improve accuracy and robustness, especially with low overlap data.

Contribution

It proposes a new registration approach integrating hard and soft assignments into an ICP variant for better handling partial overlaps.

Findings

Achieves superior accuracy over previous methods.

Demonstrates robustness with low-overlap point clouds.

Validated on public datasets with positive results.

Abstract

For the registration of partially overlapping point clouds, this paper proposes an effective approach based on both the hard and soft assignments. Given two initially posed clouds, it firstly establishes the forward correspondence for each point in the data shape and calculates the value of binary variable, which can indicate whether this point correspondence is located in the overlapping areas or not. Then, it establishes the bilateral correspondence and computes bidirectional distances for each point in the overlapping areas. Based on the ratio of bidirectional distances, the exponential function is selected and utilized to calculate the probability value, which can indicate the reliability of the point correspondence. Subsequently, both the values of hard and soft assignments are embedded into the proposed objective function for registration of partially overlapping point clouds and a novel variant of ICP algorithm is proposed to obtain the optimal rigid transformation. The proposed approach can achieve good registration of point clouds, even when their overlap percentage is low. Experimental results tested on public data sets illustrate its superiority over previous approaches on accuracy and robustness.

Figures

Tables

Table 1. Table 1: Efficiency comparison of BdICP and the proposed approach for different overlapping point clouds.

No. of scan pair	Overlap	Turb. 1		Turb. 2		Turb. 3
		BdICP	Ours	BdICP	Ours	BdICP	Ours
1	47.11%	1.9889	0.4696	1.6356	0.5643	6.4039	2.6563
2	54.92%	2.0694	0.6242	1.4216	0.8189	6.7972	2.9957
3	63.19%	2.3623	1.5701	2.7494	1.1624	4.4361	2.1182
4	72.63%	4.0802	4.1701	5.7047	3.0769	4.1523	2.0203
5	81.05%	1.7954	0.8091	1.3961	0.5694	2.4287	2.2407

Table 2. Table 2: Performance comparison of all competed approaches tested on five datasets.

Dataset	Overlap	$N_m$	$d$	$N_d$	Terms	FTrICP	wICP	CtICP	Ours
Buddha	33.29%	12147	0.0463	15233	${\epsilon }_{𝐑}$	1.0087	1.2292	0.4599	0.2015
					${\epsilon }_{\overrightarrow{t}}$	0.2560	0.3470	1.0240	0.3394
					$T\left(s\right)$	5.5287	8.4895	4.8590	0.3312
Dragon	48.36%	4418	0.0627	2976	${\epsilon }_{𝐑}$	0.2380	0.1864	0.7242	0.1042
					${\epsilon }_{\overrightarrow{t}}$	0.2752	0.2362	1.2914	0.2423
					$T\left(s\right)$	0.3228	0.4017	0.8719	0.0630
Bunny	54.51%	8051	0.0608	7205	${\epsilon }_{𝐑}$	0.2013	0.2356	0.7826	0.0765
					${\epsilon }_{\overrightarrow{t}}$	0.1714	0.2076	1.0398	0.0666
					$T\left(s\right)$	1.1260	1.4157	2.4135	0.1861
Armadillo	61.93%	3846	0.0632	5558	${\epsilon }_{𝐑}$	0.1514	0.1534	0.4626	0.0713
					${\epsilon }_{\overrightarrow{t}}$	0.0670	0.0644	0.7917	0.0685
					$T\left(s\right)$	0.4484	0.5147	2.1055	0.0871
Angel	72.26%	12999	0.5144	13540	${\epsilon }_{𝐑}$	0.0122	0.0112	0.0014	0.0014
					${\epsilon }_{\overrightarrow{t}}$	0.3873	0.3843	0.0382	0.0307
					$T\left(s\right)$	1.8085	1.7264	2.5578	1.5934

Equations

$$\begin{array}[]{l}\mathop{\min}\limits_{{\bf{R}},\vec{t},c\left(i\right)\in\left\{{1,\ldots,{N_{m}}}\right\}}\sum\limits_{i=1}^{{N_{d}}}{\left\|{{\bf{R}}{{\vec{d}}_{i}}+\vec{t}-{{\vec{m}}_{c\left(i\right)}}}\right\|}_{2}^{2}\\ \quad s.t.{{\bf{R}}^{T}}{\bf{R}}={{\bf{I}}_{3\times 3}},\det\left({\bf{R}}\right)=1\end{array}$$

$$c_k\left(i\right)=\underset{j\in \left\{1,2,\dots ,N_m\right\}}{\mathrm{argmin}}{\Vert {\mathbf{R}}_{k-1}{\vec{d}}_i+{\vec{t}}_{k-1}-{\vec{m}}_j\Vert }_2(i=1,2,\dots )$$

$$\left({\mathbf{R}}_k,{\vec{t}}_k\right)=\underset{\mathbf{R},\vec{t}}{\mathrm{argmin}}\left(\sum _{i=1}^{N_d}{\Vert \mathbf{R}{\vec{d}}_i+\vec{t}-{\vec{m}}_{c_k\left(i\right)}\Vert }_2^2\right)$$

$$\underset{\mathbf{R},\vec{t}}{\min }\sum _{i=1}^{N_d}{\Vert \mathbf{R}{\vec{d}}_i+\vec{t}-{\vec{m}}_{c\left(i\right)}\Vert }_2^2\cdot {\omega }_i\cdot p_i$$

$$c_k\left(i\right)=\underset{j\in \left\{1,2,\dots ,N_m\right\}}{\mathrm{argmin}}{\Vert {\mathbf{R}}_{k-1}{\vec{d}}_i+{\vec{t}}_{k-1}-{\vec{m}}_j\Vert }_2$$

$$l=\underset{j\in \{1,...,N_d\}}{\mathrm{argmin}}{\Vert {\vec{m}}_{c(i)}-({\mathbf{R}}_{k-1}{\vec{d}}_j+{\vec{t}}_{k-1})\Vert }_2$$

$${{w_{i}}=\left\{\begin{array}[]{l}1,{\rm{\quad if\quad}}{\vec{d}_{i}}\in{D_{\xi}}\\ 0,{\rm{\quad otherwise}}\end{array}\right.}$$

$$p_i=e^{-\gamma ({\rho }_i-1)}$$

$$({\mathbf{R}}_k,{\vec{t}}_k)=\underset{\mathbf{R},\vec{t}}{\mathrm{argmin}}(\sum _{i=1}^{N_d}{\Vert \mathbf{R}{\vec{d}}_i+\vec{t}-{\vec{m}}_{c\left(i\right)}\Vert }_2^2\cdot {\omega }_i\cdot p_i)$$

$$l=\underset{i\in \{1,...,N_d\}}{\mathrm{argmin}}{\Vert {\mathbf{R}}_{k-1}^{-1}{\vec{m}}_{c(i)}-{\mathbf{R}}_{k-1}^{-1}{\vec{t}}_{k-1}-{\vec{d}}_i\Vert }_2$$

$${\begin{array}[]{l}\mathop{\min}\limits_{\xi,{\bf{R}},\vec{t}}\left({\frac{1}{{\left|{{D_{\xi}}}\right|{\xi^{1+\lambda}}}}\sum\limits_{{{\vec{d}}_{i}}\in{D_{\xi}}}{\left\|{{\bf{R}}{{\vec{d}}_{i}}+\vec{t}-{{\vec{m}}_{c(i)}}}\right\|_{2}^{2}}}\right)\\ \quad s.t.{\rm{}}{{\bf{R}}^{T}}{\bf{R}}={{\bf{I}}_{3}},{\rm{det}}\left({\bf{R}}\right)=1\\ \quad{\rm{}}\xi\in\left[{{\xi_{\min}},1}\right],{\rm{}}{{\rm{D}}_{\xi}}\subseteq D,{\rm{}}\left|{{D_{\xi}}}\right|=\xi\left|D\right|\end{array}}$$

$$\psi (D_{\xi })=\frac{1}{\left|D_{\xi }\right|{\xi }^{1+\lambda }}\sum _{{\vec{d}}_i\in D_{\xi }}{\Vert {\mathbf{R}}_{k-1}{\vec{d}}_i+{\vec{t}}_{k-1}-{\vec{m}}_{c(i)}\Vert }_2^2$$

$${\rho }_i=\frac{{\Vert {\vec{d}}_i-{\vec{m}}_{c(i)}\Vert }_2}{{\Vert {\vec{m}}_{c(i)}-{\vec{d}}_l\Vert }_2}$$

$${\rho }_i=\frac{{\Vert {\vec{d}}_i-{\vec{m}}_{c(i)}\Vert }_2+\delta }{{\Vert {\vec{m}}_{c(i)}-{\vec{d}}_l\Vert }_2+\delta }$$

$$q_i={\omega }_i\cdot p_i$$

$$\begin{array}[]{l}J({\bf{R}},\vec{t})=\sum\limits_{i=1}^{{N_{d}}}{[{{({\bf{R}}{{\vec{d}}_{i}}-{{\vec{m}}_{c(i)}})}^{T}}({\bf{R}}{{\vec{d}}_{i}}-{{\vec{m}}_{c(i)}}){q_{i}}}\\ +2{{\vec{t}}^{T}}({\bf{R}}{{\vec{d}}_{i}}-{{\vec{m}}_{c(i)}}){q_{i}}+{{\vec{t}}^{\rm{T}}}\cdot{\vec{t}}\cdot{q_{i}}]\\ \end{array}$$

$$\frac{\partial J}{\partial \vec{t}}=\sum _{i=1}^{N_d}\left[2(\mathbf{R}{\vec{d}}_i-{\vec{m}}_{c(i)})q_i+2\vec{t}\cdot q_i\right]$$

$${\vec{t}}_k=-\frac{\sum _{i=1}^{N_d}(\mathbf{R}{\vec{d}}_i-{\vec{m}}_{c(i)})q_i}{\sum _{i=0}^{N_d}{q}_i}$$

$$J(\mathbf{R})=\sum _{i=1}^{N_d}{\Vert \mathbf{R}{\vec{x}}_i-{\vec{y}}_i\Vert }_2^2\cdot q_i$$

$$\begin{array}[]{l}{\bf{R}}=\mathop{\arg\min}\limits_{\bf{R}}{\sum\limits_{i=1}^{{N_{d}}}({{\vec{x}}_{i}^{T}{{\vec{x}}_{i}}+{\vec{y}}_{i}^{T}{{\vec{y}}_{i}}-2{\vec{x}}_{i}^{T}{\bf{R}}{{\vec{y}}_{i}}}})\cdot{q_{i}}\\ \quad=\mathop{{\mathop{\arg\max}\nolimits}}\limits_{\bf{R}}\sum\limits_{i=1}^{{N_{d}}}({{\vec{x}}_{i}^{T}{\bf{R}}{{\vec{y}}_{i}}}{q_{i}})\end{array}$$

$$\mathbf{H}=\sum _{i=1}^{N_d}{\vec{y}}_i{q}_i{\vec{x}}_i^T$$

$$[\mathbf{U},\Lambda ,\mathbf{V}]=\mathrm{SVD}(\mathbf{H})$$

$$\mathbf{R}=\mathbf{V}\mathbf{D}{\mathbf{U}}^T$$

$$\xi =0.95\times \left(0.95N-2n\right)/\left(0.95N-n\right)$$

Full text

††journal: ISPRS Journal of Photogrammetry and Remote Sensing

1 Introduction

Aiming at building up correspondences and calculating the optimal transformation between two or more point clouds, point cloud registration is of great importance in 3D model reconstruction [1], robot mapping [2, 3], object recognition [4, 5] and etc. As the most basic and essential perspective of multi-view registration, pair-wise registration has been devoted to by a large number of researchers. Among the great deal of methods, the most original and popular one is the iterative closest point (ICP) algorithm proposed by P. J. Besl and N. D. McKay [6], which iteratively builds up correspondences and calculates the rigid transformation by minimizing the Euclidean distance based residual error. However, this approach is time-consuming and only apply to two completely overlapping clouds.

In the past decade, a great deal of researchers have made efforts to improve the performance of the traditional ICP algorithm. To improve the low efficiency of ICP algorithm, Fitzgibbon [7] employed the Levenberg—Marquardt algorithm to speed it up; Jost [8] proposed a coarse-to-fine multi-resolution technique with the neighbor search algorithm to improve its efficiency. For the partially overlapping problem neglected in original ICP, a straightforward solution is to exclude point correspondences with distances greater than the user-specified threshold [9]. After that, Godin [10] proposed weighting of point correspondences to assign lower weights to correspondences with greater point-to-point distances. But this approach is ineffective. Further, Chetverikov et al.[11] proposed the trimmed ICP (TrICP) algorithm, in which an overlapping parameter was introduced into the residual error. This approach can trim the outliers automatically and obtain accurate registration results for partially overlapping clouds, but it is time-consuming. To solve this problem, Phillips et al.[12] proposed the fractional TrICP (FTrICP) algorithm to simultaneously compute the overlapping parameter and transformation in a much faster speed.

Although the aforementioned methods get substantial improvement in both robustness and efficiency, they still suffer either poor accuracy or local minimum problem. On one hand, to obtain more precise transformation, some probabilistic methods [13, 14, 15, 16] were proposed for partial registration of point clouds. But the huge computational resources they require pose a great challenge to most application areas. On the other hand, to avoid the local minimum problem existing in many algorithms, Lee et al. [17] introduced invariant features into ICP to decrease the chance being trapped into local minimum. Besides, the genetic algorithm [18, 19] and particle filter [20] were also utilized to determine the optimal transformation.

In order to achieve more robust and accurate registration, some original methods based on the traditional ICP were proposed. For example, Zhu [21] proposed the concept of bidirectional distances to solve the non-rigid registration of absolutely overlapping point clouds. Afterwards, backward distance measurement was introduced into affine ICP method by Ma et al.[22]. Although the bilateral correspondence based algorithms are more accurate than most other methods, their efficiency is unsatisfactory. Recently, correntropy [23, 24] was introduced into ICP algorithm and can be solved by maximizing the correntropy between two point clouds. Although it can deal with rigid registration with noises and outliers, its robustness should be further improved for registration of cloud pair with low overlap percentage.

However, all the above methods can only be solved under an initial transformation that is close to the ground truth. Therefore, a group of researchers contribute to feature extraction to provide initial transformation for pair-wise registration. These algorithms are mainly classified into two categories, i.e., global feature based and local feature based algorithms. The objective of global feature based algorithms is to construct a set of features to encode the geometric properties of the entire 3D object, such as geometric 3D moments [25], shape distribution [26] and spherical harmonics [27]. On the other hand, the objective of local feature based algorithms is to define a set of features to encode the characteristics of the local neighborhood of the feature points. For example, a set of point feature histograms was introduced to provide good starting points for ICP algorithm in [28] and [29]. Besides, Guo [30] designed a 3D object recognition algorithm based on his proposed local reference frame and rotational projection statistics. However, these algorithms usually suffer from either low descriptiveness or weak robustness.

Considering the aforementioned problems, we therefore propose an extended version of [31] to achieve more robust and efficient registration. Compared with the original version, this approach introduces a variable of hard assignment into the objective function, in which each point in the data shape is assigned a binary value to indicate its probability as an inlier. Based on this step, the variable of soft assignment is further explored to illustrate the reliability of the established point correspondence. To compare the proposed approach with that without hard assignment, some experiments were performed on some real datasets and the experimental results demonstrate the advantages over the former method in efficiency, accuracy and robustness.

The rest of this paper is organized as the following: Section 2 briefly reviews the ICP algorithm. The motivation behind our proposed approach is explained in section 3. Then in section 4, the proposed approach is illustrated in detail. After that is section 5, in which the superiority of the proposed method is further validated and some convincing experiments are demonstrated to compare the proposed algorithm with other related ones. At last, some conclusions are drawn in section 6.

2 The ICP algorithm

Given two point clouds, dnoted as the data shape $D=\{{\vec{d}_{i}}\}_{i=1}^{{N_{d}}}({N_{d}}\in{\mathbb{N}^{3}})$ and the model shape $M=\{{\vec{m}_{j}}\}_{j=1}^{{N_{m}}}({N_{m}}\in{\mathbb{N}^{3}})$, the goal of registration is to find the optimal rigid transformation $({\bf{R}},\vec{t})$, with which $D$ can be in the best alignment with $M$. Accordingly, it can be formulated as the following least square (LS) problem:

[TABLE]

where ${\bf{R}}\in\mathbb{R}^{3\times 3}$ is the rotation matrix, $\vec{t}\in{\mathbb{R}^{3}}$ indicates the translation vector, ${\vec{m}_{c(i)}}$ represents the correspondence of ${\vec{d}_{i}}$ in the model shape, and ${\left\|\right\|}_{2}$ denotes $L_{2}$ norm. Eq. (1) can be solved by the ICP algorithm [6], which achieves the rigid registration by iterations. Given an initial transformation $({\bf{R}_{0}},{\vec{t}_{0}})$, two steps are included in each iteration:

According to $\left({k-1}\right)th$ transformation, establish correspondence for each point ${\vec{d}_{i}}$ in the data shape:

[TABLE] 2. 2.

Based on the established correspondences, calculated the transformation $\left({{\bf{R}},\vec{t}}\right)$ by minimizing the following function:

[TABLE]

Finally, the optimal transformation can be obtained by repeating steps (1)-(2) until some stop conditions are satisfied. Although the ICP algorithm has good performances for the rigid registration, it cannot achieve the registration of partially overlapping point clouds.

3 Motivation

For point cloud registration, two issues should be carefully considered in practical applications: 1) The point cloud to be registered always contains non-overlapping areas, which should be confirmed and discarded during registration; 2) Even in the overlapping areas, there seldom exits real point correspondences due to the limited sensor resolution or noise.

Given two partially overlapping point clouds, bilateral correspondences of one point in the data shape can be established by applying the search method of nearest neighbor twice. As shown in Fig. 1, for each point ${\vec{d}}_{i}$ in the data shape $D$, its nearest neighbour ${\vec{m}}_{c}(i)$ with the forward distance can be searched from the model shape $M$. Then, for each of these nearest neighbours in $M$, its correspondence ${\vec{d}}_{l}$ with the backward distance can also be searched from $D$. Consequently, two interesting phenomena can be reasonably discovered.

Observation 1: The forward distance of the point in overlapping areas is smaller than that of the point in non-overlapping areas: For one point of the data shape, if it is located in the overlapping areas, its correspondence can be found in the overlapping area of the model shape. Therefore, the point pair of an established correspondence is very close to each other and the forward distance is always small. However, for one point located in the non-overlapping areas, it has no geometrically consistent correspondence. Hence, the point pair of an established correspondence is much far from each other and its forward distance is always larger than that of the point located in the overlapping area. Accordingly, the value of forward distance can be viewed as the index to indicate whether one point is located in the overlapping areas or not.

Observation 2: The difference between forward and backward distances can denote the reliability of the established point correspondence. For the real point correspondence, the backward distance should be equal to its forward distance. In practice, the two points of one correspondence do not have the same position on the object due to the limited senor resolution and noise. For one point in the data shape, the backward distance is usually shorter than the forward one. And there is only a small number of points whose forward distances are equal to their backward distances. As shown in Fig. 1(b) and (c), the less the difference between forward and backward distances is, the more likely this point pair of an established correspondence truly matches. In other words, with the increase of the difference, the reliability of point correspondence decreases. Therefore, the relationship of bidirectional distances can be viewed as an index to indicate whether the established point correspondence is reliable or not. If bidirectional distances of one point are almost equal, then its established correspondence can be viewed as the real correspondence with a higher probability. Otherwise the probability should accordingly decrease.

Actually, the points located in the non-overlapping areas are outliers, which should not be considered in registration. Besides, as the reliability of each point correspondence varies, different emphasis should be placed on individual point correspondence. Accordingly, a novel approach can be motivated to achieve the registration of partially overlapping point clouds.

4 The proposed approach

In this section, an effective approach is proposed for registration of partially overlapping point clouds.

4.1 ICP algorithm based on the hard and soft assignments

Based on the two observations, two variables can be introduced into point cloud registration: the variable for the hard assignment and the variable for the soft assignment. The former variable is a binary parameter utilized to indicate whether one point is an inlier or outlier, while the latter one is set to be a continuous variable between 0 and 1, which can denote the reliability of the established point correspondence. By introducing the former variable, outliers can be rejected during registration. Moreover, the latter variable is introduced to place more emphasis on these established point correspondences with high reliabilities.

Accordingly, a new objective function based on both hard and soft assignments can be proposed for registration of partially overlapping point clouds:

[TABLE]

in which ${\omega_{i}}$ is the probability denoting ${\vec{d}_{i}}$ as an inlier and ${p_{i}}$ is the reliability of the established correspondence ${\vec{d}_{i}}$ and ${\vec{m}_{c\left(i\right)}}$. The binary probability ${\omega_{i}}$ is easy to determine since there is only two cases for one point, i.e. inlier or outlier. However, the definition of ${p_{i}}$ is much more complex considering its property as a real parameter. To take full advantage of observation 2, we turn to an exponential function based probability model to design accurate reliability of an established correspondence.

Noting that Eq. (4) is an extension of Eq. (1), a similar solving approach to the traditional ICP algorithm can be adopted by iterating the following three steps:

Based on $\left({k-1}\right)th$ transformation, establish bilateral correspondences for each point ${\vec{d}_{i}}$ in $D$:

[TABLE]

[TABLE] 2. 2.

Calculate the variables of hard assignment and soft assignment, respectively:

[TABLE]

[TABLE]

in which ${D_{\xi}}$ is a subset of the data shape representing the overlapping areas in $D$, ${\rho_{i}}$ is the ratio of bidirectional distances and $\gamma$ is a preset parameter. 3. 3.

Update the rigid transformation by minimizing the objective function:

[TABLE]

Under the given initial transformation, the optimal solution can be obtained by iteratively performing the above three steps until convergence criteria are satisfied.

4.2 Establishing the bilateral correspondences

In this paper, the $k-d$ tree based nearest neighbor searching method is adopted to find correspondences. For the transformed point $({{\bf{R}}_{k-1}}{\vec{d}_{i}}+{\vec{t}_{k-1}})$ in Eq. (5), its nearest neighbor can be searched from the model shape $M$ in each iteration. Therefore, the $k-d$ tree for $M$ can be built before starting iteration and it only requires to be built once. However, the $k-d$ tree for the transformed data shape ${T_{k-1}}(D)$ requires be continuously rebuilt in each iteration since the nearest neighbor for ${\vec{m}_{c(i)}}$ is searched from the transformed data shape each time, resulting in low efficiency of the whole approach. To improve efficiency, a modified version of Eq. (6) is illustrated as:

[TABLE]

Obviously, Eq. (6) and (10) can achieve the same result. In this way, the search process can be accomplished through the following two steps: firstly, the $k-d$ tree based search method is adopted to find the nearest neighbor ${\vec{m}_{c(i)}}$ for each point in $D$; subsequently, for all the transformed nearest neighbors $({\bf{R}}_{k-1}^{-1}{\vec{m}_{c(i)}}-{\bf{R}}_{k-1}^{-1}{\vec{t}_{k-1}})$, their corresponding nearest neighbors are searched from the data shape. Therefore, the $k-d$ tree for $D$ can be built only once before starting iteration.

In practical implementation, there is no need establishing bilateral correspondences for every point in $D$. Since the hard assignment has excluded most outliers from registration, we can only search bilateral correspondences for each inlier in the data shape, as shown in Fig. 1(b) and (c). Therefore, computational efficiency can be improved significantly.

4.3 Calculating the variables of hard and soft assignments

The aim of introducing hard assignment into registration is to assign a binary value for each point in $D$ to indicate its probability as an inlier. Considering that several parameters need determining in Eq. (4), we thus turn to a simpler but similarly accurate method—the TrICP algorithm—for solution, in which the objective function is formulated as:

[TABLE]

where $\lambda$ is a preset parameter, $\left|\cdot\right|$ denotes the cardinality of a set and $\xi$ is the overlapping percentage. Accordingly, those points that minimize the above function are viewed as inliers while others are outliers.

On the basis of [12], Eq. (11) can be solved in a sequence processing manner by sorting the point correspondences $\left\{{{{\vec{m}}_{c(i)}},({{\bf{R}}_{k-1}}{{\vec{d}}_{i}}+{{\vec{t}}_{k-1}})}\right\}_{c(i)=1}^{{N_{m}}}$ according to their distances in ascending order and calculating

[TABLE]

At each iteration, a pair of sorted points is added to compute the corresponding value of $\psi({D_{\xi}})$. By traversing all sorted point correspondences, it is easy to obtain the minimum value $\psi({D_{\xi}})$ that corresponds to the optimal overlapping areas ${D_{\xi}}$. Subsequently, points belonging to ${D_{\xi}}$ are assigned the probability $1$.

Considering the property of the difference between forward and backward distances for a point in the data shape, we therefore adopt Eq. (8) as the variable of soft assignment. We originally design the ${\rho_{i}}$ as the ratio of bidirectional distances, which is formulated as:

[TABLE]

in which ${\vec{d}_{l}}$ is the backward correspondence of ${\vec{m}_{c(i)}}$, ${\vec{d}_{i}}=({{\bf{R}}_{k-1}}{\vec{d}_{i}}+{\vec{t}_{k-1}})$ and ${\vec{d}_{l}}=({{\bf{R}}_{k-1}}{\vec{d}_{l}}+{\vec{t}_{k-1}})$.

Usually, the ratio of bidirectional distances is expected to meet the condition ${\rho_{i}}\in[1,+\infty)$. Sometimes, the point correspondence established by Eq. (5) or Eq. (6) may be coincident due to the round-off effect. In that case, ${\rho_{i}}$ will be illegal. Therefore, Eq. (13) should be modified as follows:

[TABLE]

where $\delta$ is a small positive parameter.

In practical applications, the effects of hard and soft assignments rely on parameters such as $\lambda$ and $\gamma$. For accurate registration, we want to select as many as inliers, so some outliers may be included into registration. To eliminate the disturbance, the parameter $\gamma$ can be adjusted to a smaller value to give hard penalty to these selected outliers.

4.4 Computation of the rigid transformation

To compute the rigid transformation $({\bf{R}},\vec{t})$, the following derivation presents the closed-form solution for this question.

Substitute ${q_{i}}$ for the product of two variables of hard and soft assignments:

[TABLE]

Extend the quadratic-formed objective function depicted in Eq. (4) as following:

[TABLE]

Take the derivative of $J$ with respective to $\vec{t}$, the following result can be acquired:

[TABLE]

Let $\frac{{\partial J}}{{\partial{\vec{t}}}}=0$, the translate vector can be obtained as:

[TABLE]

In Eq. (16), $\vec{t}$ can be replaced by Eq. (18) and the objective function can be simplified as:

[TABLE]

where ${\vec{x}_{i}}={\vec{d}_{i}}-{{\sum\limits_{i=1}^{{N_{d}}}{{q_{i}}}{{\vec{d}}_{i}}}\mathord{\left/{\vphantom{{\sum\limits_{i=1}^{{N_{d}}}{{q_{i}}}{{\vec{d}}_{i}}}{\sum\limits_{i=1}^{{N_{d}}}{{q_{i}}}}}}\right.\kern-1.2pt}{\sum\limits_{i=1}^{{N_{d}}}{{q_{i}}}}}$ and ${\vec{y}_{i}}={\vec{m}_{c(i)}}-{{\sum\limits_{i=1}^{{N_{d}}}{{q_{i}}}{{\vec{m}}_{c(i)}}}\mathord{\left/{\vphantom{{\sum\limits_{i=1}^{{N_{d}}}{{q_{i}}}{{\vec{m}}_{c(i)}}}{\sum\limits_{i=1}^{{N_{d}}}{{q_{i}}}}}}\right.\kern-1.2pt}{\sum\limits_{i=1}^{{N_{d}}}{{q_{i}}}}}$. Subsequently, the rotation matrix can be calculated by minimizing the function $J(\bf{R})$, which can be expanded as:

[TABLE]

Taking advantage of the optimization for Eq. (20) solved by Myronenko [32], the conclusion can be directly presented as:

1. Compute the $3\times 3$ matrix $\bf{H}$ and its singular value decomposition (SVD) results:

[TABLE]

[TABLE]

2. Calculate the rotation matrix:

[TABLE]

where ${\bf{D}}=diag(1,1,\det({\bf{V}}{{\bf{U}}^{T}}))$.

4.5 Algorithm implementation

According to the above analysis and derivation, the implementation of the proposed approach can be summarized in Algorithm 1.

Given the initial parameters $({{\bf{R}}_{0}},{\vec{t}_{0}})$, the proposed approach can achieve registration of partially overlapping point clouds with good performances.

5 Experimental results

In this section, we firstly validated the necessity of the hard and soft assignments combined method for registration. Through a series of experimental results, the superiority of introducing hard assignment into registration is self-evident. Then in the second part, the proposed approach was compared with several methods both for data simulation and for experiment comparison. Finally, the proposed method was applied to some large-scale environmental point clouds. All of the following experiments adopt the $k-d$ tree based nearest neighbor searching method on one same computer.

5.1 Validation

A significant advantage of introducing soft assignment into registration is the high accuracy it brings about. However, besides the low efficiency due to bilateral correspondences for all points in the data shape, it can not avoid convergence of local minimum neither. By introducing hard assignment to exclude outliers, this phenomenon can be alleviated remarkably. Fig. 2 displays the registration results without and with the hard assignment strategy. For a pair of point clouds from the Stanford Bunny [33], the algorithm using soft assignment converges to local minimum and results in failing registration, while our proposed method can obtain precise results, as shown in Fig. 2 (a) and (b), respectively.

To further demonstrate the necessity of introducing hard assignment, the following experiments were designed to compare the performances without and with the hard assignment. In this experiment, five pairs of point clouds scanned from different viewpoints of one same object were randomly chosen. For each scan pair, the ground truth of transformation $\left({{{\bf{R}}_{g}},{{\vec{t}}_{g}}}\right)$ was imposed by different perturbed transformations $\left({{{\bf{R}}_{r}},{{\vec{t}}_{r}}}\right)$, which were randomly drawn from the uniformation distribution $U\left[{-5\deg,5\deg}\right]$, $U\left[{-5{\mathop{\rm d}\nolimits},5{\mathop{\rm d}\nolimits}}\right]$; $U\left[{-10\deg,10\deg}\right]$, $U\left[{-5{\mathop{\rm d}\nolimits},5{\mathop{\rm d}\nolimits}}\right]$; and $U\left[{-10\deg,10\deg}\right]$, $U\left[{-7{\mathop{\rm d}\nolimits},7{\mathop{\rm d}\nolimits}}\right]$ respectively, in which $d$ is the average point resolution of model shape. We respectively denote these turbulence as Turb. 1, Turb. 2 and Turb. 3. Taking the rigid transformation $\left({{{\bf{R}}_{r}}\cdot{{\bf{R}}_{g}},{{\vec{t}}_{r}}+{{\vec{t}}_{g}}}\right)$ as initial transformation, the estimated transformation $\left({{{\bf{R}}_{m}},{{\vec{t}}_{m}}}\right)$ can be obtained accordingly. For comparing accuracy of both approaches, relative errors are defined as: ${\varepsilon_{\bf{R}}}={\left\|{{{\bf{R}}_{m}}-{{\bf{R}}_{g}}}\right\|_{F}}$ and ${\varepsilon_{\vec{t}}}=\frac{{{{\left\|{{{\vec{t}}_{m}}-{{\vec{t}}_{g}}}\right\|}_{2}}}}{d}$. Table 1 records the overlapping percentage of each scan pair and the runtime for the methods before and after introducing hard assignment under different turbulences. To demonstrate the results figuratively, Fig. 3 displays the registration results.

As demonstrated in Fig. 3, in each trial, the relative errors generated from the proposed method are smaller than that without using hard assignment. Although the errors of both methods will increase with the raise of turbulence, the proposed method merely gains a moderate increase while the method solely using soft assignment has an apparent augmentation. As for the soft assignment approach, the registration results for higher overlapping cloud pairs were not always better than that for lower overlapping cloud pairs under the same turbulence. This is because the approach without using hard assignment is easy to be trapped into local minimum and thus always results in failing registration. Best illustrated in low overlapping situations, the proposed algorithm could effectively avoid the local minimum problem. Besides, Table 1 demonstrates a distinct comparison about efficiency. As it shows, the efficiency of the proposed method showed an obvious superiority over merely soft assignment method, especially when the overlap percentage was low. With the increase of turbulences, efficiency of both methods will decrease, but the proposed method still outperformed the method without using hard assignment. This high efficiency owes to the hard assignment strategy, which excludes most outliers from establishing bilateral correspondences and computing the ratio of bilateral distances. Obviously, the proposed method using both hard and soft assignments is superior to that merely using soft assignment in efficiency, accuracy and robustness.

5.2 Comparison

To testify its superior performances, the proposed approach was compared with the fractional TrICP algorithm (FTrICP), ICP algorithm with weighting of pairs (wICP) [9] and the correntropy ICP algorithm (CtICP) on several real datasets. For providing solid proof about their performances, the experiments were conducted both for data simulation and for experiment comparison.

5.2.1 Data simulation

To perform the simulation comparison, a point cloud with $8171$ 3D points was selected from the Stanford Bunny. For generating the data shape, the following processes were implemented: The original shape was randomly deleted $5\%$ points and then cut down by one part, which included $n$ range points; then the remaining shape was added with Gaussian zero mean noise. For generating the model shape, the following two steps were implemented: Deleting $5\%$ points in the original shape and applying a randomly generated transformation $\left({{{\bf{R}}_{r}},{{\vec{t}}_{r}}}\right)$ to the remaining shape; then cutting down the transformed shape by other $n$ range points. Accordingly, the overlap percentage of the two generated shapes is calculated as:

[TABLE]

in which $N$ denotes the number of points of the original shape. To simulate the registration results under different overlap percentages, the value $n$ can be accordingly changed. Thus we can obtain a set of registration results for all competed approaches.

For comparison, relative errors are defined as ${\varepsilon_{\bf{R}}}={\left\|{{{\bf{R}}_{m}}-{{\bf{R}}_{r}}}\right\|_{F}}$ and ${\varepsilon_{\vec{t}}}={\left\|{{{\vec{t}}_{m}}-{{\vec{t}}_{r}}}\right\|_{2}}$, where $\left({{{\bf{R}}_{m}},{{\vec{t}}_{m}}}\right)$ represents the estimated transformation, and ${\left\|{}\right\|_{F}}$ denotes Frobenius norm. Denote $d$ as the average point resolution of model shape, the registration process can be viewed as successful if and only if (${\varepsilon_{\bf{R}}}\leq 0.01$ and ${\varepsilon_{\vec{t}}}\leq d$). All approaches were tested on each cloud pairs with varied overlap percentages for $100$ Monte Carlo (MC) trials. To view it more intuitively, Fig. 4 displays registration performances for all competed approaches.

As shown in Fig. 4(a), both wICP and CtICP are satisfactory for pretty highly overlapping clouds. This is because the correspondence distance threshold specified by user is mainly determined empirically and many outliers may be included into registration. The FTrICP algorithm can achieve more superior performance compared with the former two methods. But its robustness decreases quickly with the decrease of overlap percentage. This is a common phenomenon since many methods neglect the reliability of established correspondences. On the contrary, the proposed approach can still obtain precise registration results even though the overlap percentage is pretty small. Moreover, its robustness decreases quite slowly with the increase of outliers. Benefited from the hard assignment, this method can gain a slight raise on registration precision when the overlap percentage is small enough. Therefore, the proposed approach achieves the most robust registration among all the competed approaches. In terms of efficiency, Fig. 4(b) gives the best illustration. When the overlap percentage is pretty high, the proposed approach has the lowest efficiency since it needs to compute bilateral distances for most points. Once the overlap percentage is lower than $80\%$, its superiority over other approaches becomes distinct and even has the highest efficiency among all the competed approaches. This benefits from the hard assignment strategy, which has excluded most outliers from registration.

5.2.2 Experiment comparison

In this section, the proposed approach was compared with the aforementioned methods on several real datasets: Happy Buddha, Dragon, Bunny, Armadillo and Angel. Detailed information of each cloud pair was depicted in Table 2. In addition to point clouds information, the ground truth of rigid transformation $\left({{{\bf{R}}_{g}},{{\vec{t}}_{g}}}\right)$ for each pair of cloud is also known.

Before registration, a perturbed rigid transformation $\left({{{\bf{R}}_{r}},{{\vec{t}}_{r}}}\right)$ was imposed to the ground truth, in which the Euler angle and translate component were randomly drawn from the uniform distribution $U\left[{-5\deg,5\deg}\right]$ and $U\left[{-{\mathop{\rm d}\nolimits},{\mathop{\rm d}\nolimits}}\right]$ respectively. Here, the symbol $d$ denotes the average point resolution of the model shape. Taking the rigid transformation $\left({{{\bf{R}}_{r}}\cdot{{\bf{R}}_{g}},{{\vec{t}}_{r}}+{{\vec{t}}_{g}}}\right)$ as initial parameter for all competed approaches, the estimated transformation $\left({{{\bf{R}}_{m}},{{\vec{t}}_{m}}}\right)$ can be obtained accordingly. For comparing accuracy of each approach, relative errors are defined as: ${\varepsilon_{\bf{R}}}={\left\|{{{\bf{R}}_{m}}-{{\bf{R}}_{g}}}\right\|_{F}}$ and ${\varepsilon_{\vec{t}}}=\frac{{{{\left\|{{{\vec{t}}_{m}}-{{\vec{t}}_{g}}}\right\|}_{2}}}}{d}$. All competed approaches were tested on each dataset with $20$ MC trials to eliminate random interference. To evaluate the performance of each approach, Table 2 demonstrates the relative error of rotation matrix, translation component and average runtime, respectively. To view this in a more intuitive way, Fig. 5 displays the registration results for all competed approaches.

As shown in Table 2 and Fig. 5, all approaches can achieve precise registration on Angel, which has high overlap percentage. Since it neglects outliers and relies highly on initial transformation, the CtICP algorithm obtained unsatisfactory results on the other four datasets. While FTrICP and wICP achieved more precise and efficient registration than CtICP on most datasets, they still fall behind our proposed approach by a large margin. As for the wICP algorithm, it neglects the nonoverlapping areas problem and searches a correspondence for each point in the data shape, which includes many outliers into registration. Although the FTrICP algorithm takes outliers into consideration, it deems that all the established correspondences are correct, thus may fail when registration are interfered by many outliers. The relative errors of ${\bf{R}}$ and ${\vec{t}}$ increase dramatically for the FTrICP, wICP and CtICP with the decrease of overlap percentage, but ours has a moderate rise profited from the combination of hard and soft assignments. Besides, the high efficiency of the proposed method resulted from the hard assignment also put it in an advantageous status. Based on these analysis, it is fair to draw that the proposed method is the most robust and efficient one among all these competed approaches.

5.3 Application

To confirm its application to environmental scenes, the proposed approach was testified in the following large-scale environmental point clouds, which were recorded in the Institute of Robotics and Mechatronics at the German Aerospace Center (DLR) [34]. For conducting this experiment, the given pose was imposed by a perturbed transformation $\left({{{\bf{R}}_{r}},{{\vec{t}}_{r}}}\right)$, and then used as the initial transformation. Note that the given pose was obtained via some methods and are available online. Under such deployment, the registration results using the proposed approach can be demonstrated in Fig. 6.

As illustrated in Fig. 6, the proposed approach is appliable to environmental point clouds and is able to achieve pretty precise registration, even if the point cloud contains tens of millions points. This is easy to interpret because of the hard and soft assignments combined method. On one hand, the hard assignment exclude most outliers from registration. Therefore, the reconstruction of the environmental scene is mainly based on inliers, which are essential for registration. On the other hand, the soft assignment strategy assigns a value for each pair of correspondence to denote its reliability of the established correspondence. This method takes the probability into consideration and overcomes the disadvantage of merely optimizing the least square function, thus can achieve more precise registration.

6 Conclusion

This paper proposes a novel approach for registration of partially overlapping point clouds. The novelty of this approach is the introduction of both the hard and soft assignments, which can be represented as a binary variable and continuous variable, respectively. By consideration of the forward distance, each point in the data shape can be assigned a binary value, which can denote whether it is an inlier or outlier. Meanwhile, based on the ratio of bidirectional distances, each established point correspondence can be assigned a continuous probability, which can be viewed as an index to indicate the reliability of point correspondence. By integrating both the hard and soft assignments, a new objective function can be designed for registration of partially overlapping point clouds and a variant of the ICP algorithm is proposed to calculate the optimal rigid transformation. The proposed approach has been tested on public datasets and experimental results illustrate the proposed approach can achieve the registration of partially overlapping clouds with good robustness and accuracy.

Acknowledgements

This work is supported by the National Natural Science Foundation of China under Grant Nos. 61573273, 61573280 and 61503300.