Effective scaling registration approach by imposing the emphasis on the scale factor

TL;DR

This paper introduces a novel scaling registration method that emphasizes the scale factor to avoid ill-posed problems, using a new objective function and an iterative algorithm, demonstrating improved efficiency and robustness in practical applications.

Contribution

It proposes a new objective function for scaling registration that emphasizes the scale factor, along with an iterative algorithm to improve accuracy and robustness.

Findings

Outperforms previous methods in efficiency and robustness

Effective in aligning partially overlapping point sets

Successfully applied to merging grid maps of different resolutions

Abstract

This paper proposes an effective approach for the scaling registration of $m$-D point sets. Different from the rigid transformation, the scaling registration can not be formulated into the common least square function due to the ill-posed problem caused by the scale factor. Therefore, this paper designs a novel objective function for the scaling registration problem. The appearance of this objective function is a rational fraction, where the numerator item is the least square error and the denominator item is the square of the scale factor. By imposing the emphasis on scale factor, the ill-posed problem can be avoided in the scaling registration. Subsequently, the new objective function can be solved by the proposed scaling iterative closest point (ICP) algorithm, which can obtain the optimal scaling transformation. For the practical applications, the scaling ICP algorithm is further extended to align partially overlapping point sets. Finally, the proposed approach is tested on public data sets and applied to merging grid maps of different resolutions. Experimental results demonstrate its superiority over previous approaches on efficiency and robustness.

Tables

Table 1. Table 1: Performance comparison of competing approaches tested on different 2D shapes

		BiTrICP		TrICPBNs		TrICPBWs		sTrICP
Datasets	$s$	MSE	T(s)	MSE	T(s)	MSE	T(s)	MSE	T(s)
Apple	1.4733	1.4508	0.1893	0.8277	0.0436	0.8277	0.0466	0.8275	0.0352
Bat	0.75	0.2682	0.8372	7.0585	0.1979	196.0914	0.3584	0.2680	0.3771
Chicken	1.5	0.2668	0.3368	0.2675	0.1741	1.7866	0.1896	0.2675	0.0755

Table 2. Table 2: Performance comparison of all competing approaches tested on different range images

		BiTrICP		TrICPBNs		TrICPBWs		sTrICP
Datasets	$s$	MSE	T(s)	MSE	T(s)	MSE	T(s)	MSE	T(s)
Bunny	0.5	1.2676	3.7437	1.2747	1.1191	4.1380	2.9040	1.2698	1.6219
Dragon	0.25	0.5560	3.4324	9.0729	4.3301	18.2296	1.8909	0.5539	2.2583
HappyBuddha	1.25	0.6266	5.3875	0.6206	2.0002	0.6206	1.9986	0.6280	1.9558

Equations

$$\begin{array}[]{*{20}{c}}{\mathop{\min}\limits_{\;{\bf{R}},\vec{t},c(i)\in\{1,...,{N_{q}}\}}\sum\limits_{i=1}^{{N_{p}}}{\left\|{{\bf{R}}{{\vec{p}}_{i}}+\vec{t}-{{\vec{q}}_{c(i)}}}\right\|_{2}^{2}}}\\ {{\rm{s}}{\rm{.t}}{\rm{.}}{{\bf{R}}^{T}}{\bf{R}}={{\bf{I}}_{m}},\det({\bf{R}})=1}\\ \end{array},$$

$$c_k(i)=\underset{j\in \{1,2,...,N_q\}}{\mathrm{argmin}}{\Vert {\mathbf{R}}_{\mathbf{k}-1}{\vec{p}}_i+{\vec{t}}_{k-1}-{\vec{q}}_j\Vert }_2$$

$$({\mathbf{R}}_k,{\vec{t}}_k)=\underset{\mathbf{R},\vec{t}}{\mathrm{argmin}}\sum _{i=1}^{N_p}{\Vert \mathbf{R}{\vec{p}}_i+\vec{t}-{\vec{q}}_{c_k(i)}\Vert }_2^2$$

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

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

$$c_k(i)=\underset{j\in \{1,2,..,N_q\}}{\mathrm{argmin}}{\Vert s_{k-1}{\mathbf{R}}_{k-1}{\vec{p}}_i+{\vec{t}}_{k-1}-{\vec{q}}_j\Vert }_2\mathrm{for}i=1,2,...,N_p$$

$$(s_k,{\mathbf{R}}_k,{\vec{t}}_k)\underset{s,\mathbf{R},\vec{t}}{=\arg \min }\frac{{\sum }_{i=1}^{{}_{{}_{N_p}}}{\Vert s\mathbf{R}{\vec{p}}_i+\vec{t}-{\vec{q}}_{c_k(i)}\Vert }_2^2}{s^2}$$

$$\vec{t}=\frac{1}{N}\sum _{i=1}^{N_p}{\vec{q}}_{c_k(i)}-\frac{1}{N}\sum _{i=1}^{N_p}s\mathbf{R}{\vec{p}}_i.$$

$$F(s,\mathbf{R})=\sum _{i=1}^{{}_{{}_{N_p}}}{\Vert s\mathbf{R}({\vec{p}}_i-\frac{1}{N}{\sum }_{i=1}^{N_p}{\vec{p}}_i)-({\vec{q}}_{c_k(i)}-\frac{1}{N}{\sum }_{i=1}^{N_p}{\vec{q}}_{c_k(i)})\Vert }_2^2/s^2.$$

$$\begin{array}[]{l}F(s,{\bf{R}})=\sum\limits_{i=1}^{{N_{p}}}{\left\|{s{\bf{R}}{{\vec{d}}_{i}}-{{\vec{m}}_{i}}}\right\|_{2}^{2}}/{s^{2}}\\ \quad\quad\quad=({s^{2}}\sum\limits_{i=1}^{{N_{p}}}{\vec{d}_{i}^{T}}{{\vec{d}}_{i}}-2s\sum\limits_{i=1}^{{N_{p}}}{\vec{d}_{i}^{T}}{\bf{R}}{{\vec{m}}_{i}}+\sum\limits_{i=1}^{{N_{p}}}{\vec{m}_{i}^{T}}{{\vec{m}}_{i}})/{s^{2}}.\\ \end{array}$$

$$\mathbf{H}=\frac{1}{N_p}\sum _{i=1}^{N_p}{\vec{d}}_i^T{\vec{m}}_i$$

$$U\Lambda V=svd(\mathbf{H})$$

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

$$\frac{\partial F(s,\mathbf{R})}{\partial s}=2s^{-2}\sum _{i=1}^{N_p}{\vec{d}}_i^T\mathbf{R}{\vec{m}}_i-2s^{-3}\sum _{i=1}^{N_p}{\vec{m}}_i^T{\vec{m}}_i.$$

$$s_k=\frac{{\sum }_{i=1}^{N_p}{\vec{m}}_i^T{\vec{m}}_i}{{\sum }_{i=1}^{N_p}{\vec{d}}_i^T{\mathbf{R}}_k{\vec{m}}_i}$$

$$\vec{t}=\frac{1}{N}\sum _{i=1}^{N_p}{\vec{q}}_{c_k(i)}-\frac{1}{N}\sum _{i=1}^{N_p}{s}_k{\mathbf{R}}_k{\vec{p}}_i$$

$$\Psi (\xi ,s,\mathbf{R},\vec{t})=\frac{e(\xi ,s,\mathbf{R},\vec{t})}{s^2\cdot {\xi }^{1+\lambda }},$$

$$e(\xi ,s,\mathbf{R},\vec{t})=\frac{\sum _{{\vec{p}}_i\in P_{\xi }}{\Vert s\mathbf{R}{\vec{p}}_i+\vec{t}-{\vec{q}}_{c(i)}\Vert }_2^2}{N_p^{{}^{\prime }}}.$$

$$c_k(i)=\underset{j\in \{1,2,..,N_q\}}{\mathrm{argmin}}{\Vert s_{k-1}{\mathbf{R}}_{k-1}{\vec{p}}_i+{\vec{t}}_{k-1}-{\vec{q}}_j\Vert }_2\mathrm{for}i=1,2,...,N_p.$$

$$({\xi }_k,P_{{\xi }_k})\underset{\xi ,P_{\xi }}{=\arg \min }\Psi (\xi ,s,\mathbf{R},\vec{t}).$$

$$(s_k,{\mathbf{R}}_k,{\vec{t}}_k)\underset{s,\mathbf{R},\vec{t}}{=\arg \min }\sum _{{\vec{p}}_i\in P_{{\xi }_k}}{\Vert s\mathbf{R}{\vec{p}}_i+\vec{t}-{\vec{q}}_{c_k(i)}\Vert }_2^2/s^2.$$

$$e_k=\frac{\sum _{{\vec{p}}_i\in P_{{\xi }_{k-1}}}{\Vert {\vec{p}}_{i,k-1}-{\vec{q}}_{c_k(i)}\Vert }_2^2}{\left|P_{{\xi }_{k-1}}\right|\cdot {s_{k-1}}^2\cdot {({\xi }_{k-1})}^{1+\lambda }}.$$

$${\eta }_k=\frac{\sum _{{\vec{p}}_i\in P_{{\xi }_k}}{\Vert {\vec{p}}_{i,k-1}-{\vec{q}}_{c_k(i)}\Vert }_2^2}{\left|P_{{\xi }_k}\right|\cdot {s_{k-1}}^2\cdot {({\xi }_k)}^{1+\lambda }}$$

$${\epsilon }_k=\frac{\sum _{{\vec{p}}_i\in P_{{\xi }_k}}{\Vert {\vec{p}}_{i,k}-{\vec{q}}_{c_k(i)}\Vert }_2^2}{\left|P_{{\xi }_k}\right|\cdot {s_k}^2\cdot {({\xi }_k)}^{1+\lambda }}.$$

$$e_{k+1}=\frac{\sum _{{\vec{p}}_i\in P_{{\xi }_k}}{\Vert {\vec{p}}_{i,k}-{\vec{q}}_{c_{k+1}(i)}\Vert }_2^2}{\left|P_{{\xi }_k}\right|\cdot {s_k}^2\cdot {({\xi }_k)}^{1+\lambda }}.$$

$$e_{k+1}=\frac{\sum _{{\vec{p}}_i\in P_{{\xi }_k}}{\Vert {\vec{p}}_{i,k}-{\vec{q}}_{c_{k+1}(i)}\Vert }_2^2}{\left|P_{{\xi }_k}\right|\cdot {s_k}^2\cdot {({\xi }_k)}^{1+\lambda }}\le \frac{\sum _{{\vec{p}}_i\in P_{{\xi }_k}}{\Vert {\vec{p}}_{i,k}-{\vec{q}}_{c_k(i)}\Vert }_2^2}{\left|P_{{\xi }_k}\right|\cdot {s_k}^2\cdot {({\xi }_k)}^{1+\lambda }}={\epsilon }_k$$

Full text

††journal: Journal of Electronic Imaging

1 Introduction

Due to the wide application in 3D reconstruction [1, 2], shape recognition [3, 4], and robot mapping [5], point sets registration has attracted immense attention [6] in computer vision, pattern recognition, and robotics. The goal of registration is to establish correspondences between two point sets and recover the optimal transformation, which can well align one point set to the other. It has been the object of much attention since the seminal work presented in [7, 8].

To solve the registration problem, the most popular solution is the iterative closest point (ICP) algorithm [7]. In this approach, the rigid registration problem can be formulated into a least square (LS) function, which can be solved by establishing the point correspondences and calculating the rigid transformation iteratively. Although the ICP algorithm can achieve the rigid registration with good accuracy and efficiency, it cannot be applied to the registration of partially overlapping point sets. To trim the outliers, Chetverikov et al. [9] introduced an overlap parameter into the LS function, which can be solved by the trimmed ICP (TrICP) algorithm. By discarding outliers, this approach can obtain accurate registration results for partially overlapping point sets. Nevertheless, it is time-consuming to calculate the overlap parameter. To accelerate the registration, Phillips et al. [10] proposed the fractional TrICP (FTrICP) algorithm, which can compute the overlap parameter and the optimal transformation simultaneously. As the variants of ICP algorithm are well known to be susceptible to local convergence, the particle filter [11] or genetic algorithm [12] could be utilized to obtain the desired global minimum. Besides, effective features [13, 14, 15] can also be extracted and matched for the point sets to be registered so as to provide the initial parameters for registration approaches. These registration approaches may obtain accurate results for the rigid registration, but they may not achieve the scaling registration, which may exist in some practical applications.

To achieve the fine registration, Zha et al. [16] integrated the scale factor into a modified ICP algorithm, where extended signatures images were utilized to establish precise correspondence. Besides, Zinber et al. [17] also proposed to integrate the scale factor into the LS function, where the scale factor is directly estimated. But it can always obtain the unexpected registration results due to ill-posed problem inherited in the scaling registration. Therefore, Ying et al. [18] proposed an approach by applying the boundary constraint as a regularizer to the scale factor in the LS function and it may obtain accurate registration results. However, it is difficult to estimate the suitable range of scale factor in some cases, which may lead to the failure of registration. Subsequently, Zhu et al. [19] introduced the bidirectional distance measurement into the scaling registration and designed the corresponding LS function, which can be solved by the proposed scaling ICP algorithm. What’s more, Du et al. [20] extended the bidirectional distance measurement to deal with the scaling registration of partially overlapping point sets. Although this approach is very robust, it is somewhat time-consuming due to the establishment of bidirectional correspondences.

Meanwhile, some probabilistic approaches were also proposed for the scaling registration of partially overlapping point sets. By representing range images by Gaussian mixture models (GMMs), Myronenko and Song [21] considered the alignment of two range images as a probability density estimation problem, which can be solved by iteratively fitting the GMM centroids with EM algorithm. In [22], range images are also represented as GMMs, then the registration problem is viewed as that of aligning two Gaussian mixtures by minimizing their discrepancy. Although these probabilistic approaches are very accurate, they are always requiring huge computational resources. Recently, Du et al. [23] combined the ICP algorithm with EM algorithm to achieve the scaling registration. It also needs to estimate boundaries of the scale factor. Therefore, a more effective approach is required to be proposed for the scaling registration of point sets.

Compared to the rigid registration, the scale factor is the extra item in the scaling registration, which can increase the difficulty of this problem. As the criteria of least square (LS) error is an ill-posed problem in the scaling registration, we should pay more attention to the scale factor. By imposing the emphasis on the scale factor, we can design a new objective function, where the ill-posed problem can be transformed into the well-posed one for the scaling registration of absolutely overlapping point sets. Subsequently, the solution of this objective function will be proposed to obtain the optimal scaling transformation. For the practical applications, the proposed approach is then extended to deal with partially overlapping point sets. Besides, it will be applied to merging grid maps of different resolutions.

The remainder of this paper is organized as follows. Section 2 briefly reviews the original ICP algorithm. In Section 3, the proposed approach is presented for the scaling registration of $m$-D point sets. Following that is Section 4, in which the proposed approach is tested and evaluated on some public data sets. Finally, some conclusions are drawn in Section 5.

2 The ICP algorithm

Given two overlapping point sets in ${\mathbb{R}^{m}}$, the data shape $P=\{{\vec{p}_{i}}\}_{i=1}^{{N_{p}}}$ and the model shape $Q=\{{\vec{q}_{j}}\}_{j=1}^{{N_{q}}}({N_{p}},{N_{q}}\in{\mathbb{N}})$, the goal of rigid registration is to find the optimal transformation $({\bf{R}},\vec{t})$, with which $P$ can be in the best alignment with $Q$. Accordingly, it can be formulated as the following LS problem:

[TABLE]

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

(1) Establish new correspondence for each point ${{\vec{p}}_{i}}$ in the data shape $P$:

[TABLE]

(2) Update the rigid transformation by minimizing the following function:

[TABLE]

Although the original ICP algorithm has good performance, it cannot deal with the scaling registration problem. Besides, it does not fit for the registration of partially overlapping point sets.

3 The proposed approach

Suppose there are two overlapping point sets in ${\mathbb{R}^{m}}$, the data shape $P=\{{\vec{p}_{i}}\}_{i=1}^{{N_{p}}}$ and the model shape $Q=\{{\vec{q}_{j}}\}_{j=1}^{{N_{q}}}({N_{p}},{N_{q}}\in{\mathbb{N}})$. Accordingly, the goal of scaling registration to find the optimal scaling transformation $(s,{\bf{R}},\vec{t})$, with which $P$ can be in the best alignment with $Q$. Therefore, it seems that this registration problem can reasonably be formulated as the following LS function:

[TABLE]

where $s$ denotes the scale factor.

Although Eq. (4) is the formulation of the scaling registration, it is unreasonably adopted as the objective function of this registration problem due to one reason. As shown in Fig.1, if the scale factor $s$ tends to zero, there will appear the phenomenon that all points in data shape $P$ will be close to each other and find the same corresponding point in model shape $Q$. And when $s=0$, the objective function denoted by Eq. (4) can be minimized. But this is not the expected registration result. Therefore, Eq. (4) is an ill-posed problem, which can lead to unexpected registration results.

To address this issue, a novel objective function can be proposed by imposing the emphasis on the scale factor as follows:

[TABLE]

As Eq. (5) and Fig. 1 displayed, if the scale factor $s$ tends to zero, the item of denominator will be decreased to zero quickly, even when the item of the numerator may be gradually declined, which can avoid the ill-posed problem appeared in Eq. (4).

Here, the value of $v$ in the denominator item should be carefully chosen due to two reasons. Firstly, the item denoted by the denominator should be decreased more rapidly than or equal to that of the numerator, when the scale factor tends to zeros. Otherwise, it could not avoid the phenomenon that the scale factor tends to zero. As the item of the numerator denotes the quadratic error, the value of $v$ must be greater than or equal to 2. Secondly, it is required that Eq. (5) should contain the valid solution, which can be easily obtained. During the derivation, we find $v$=1 or $v$=2 can make it easy to obtain the solution, which is valid. Therefore, we can only choose $v$=2 in Eq. (5).

3.1 The scaling ICP algorithm

According to the original ICP algorithm, Eq. (5) can be solved by the proposed scaling ICP algorithm, which achieves the scaling registration by iterations. Given the initial transformation $(s_{0},{\bf{R}_{0}},{\vec{t}_{0}})$, two steps are included in each iteration:

(1) Establish correspondences between two $m$-D point sets:

[TABLE]

(2) Calculate the scaling transformation by minimizing the function:

[TABLE]

Step (1) can be solved by many efficient methods such as the nearest neighbor search method based on the $k$-d tree or its variants [24]. Therefore, Step (2) is the critical step.

To calculate the scaling transformation, the following lemma can be presented without proof.

Lemma 1 Given two $m$-D point sets, $D\buildrel\Delta\over{=}\{{\vec{d}_{i}}\}_{i=1}^{N}$ and $M\buildrel\Delta\over{=}\{{\vec{m}_{i}}\}_{i=1}^{N}$, then the function $F(t)=\sum\nolimits_{i=1}^{N}{\left\|{{{\vec{d}}_{i}}+\vec{t}-{{\vec{m}}_{i}}}\right\|_{2}^{2}}$ has the minimum value when $\vec{t}=\frac{1}{N}\sum\nolimits_{i=1}^{N}{{{\vec{m}}_{i}}}-\frac{1}{N}\sum\nolimits_{i=1}^{N}{{{\vec{d}}_{i}}}$.

According to this lemma, minimizing $F(s,{\bf{R}})=\sum\nolimits_{i=1}^{{}_{{}_{{N_{p}}}}}{\left\|{s{\bf{R}}{{\vec{p}}_{i}}+\vec{t}-{{\vec{q}}_{{c_{k}}(i)}}}\right\|_{2}^{2}/{s^{2}}}$ can obtain the following result:

[TABLE]

Hence, the function $F(s,{\bf{R}})$ can be transformed into the following form:

[TABLE]

Denote ${\vec{d}_{i}}={\vec{p}_{i}}-\frac{1}{{{N_{p}}}}\sum\nolimits_{i=1}^{{N_{p}}}{{{\vec{p}}_{i}}}$ and ${\vec{m}_{i}}={\vec{q}_{{c_{k}}(i)}}-\frac{1}{{{N_{p}}}}\sum\nolimits_{i=1}^{{N_{p}}}{{{\vec{q}}_{{c_{k}}(i)}}}$. Then, $F(s,{\bf{R}})$ can be simplified and expanded as:

[TABLE]

Subsequently, the scaling transformation can be calculated as follows.

(a) Rotation matrix

As shown in Eq. (10), only the middle item of $F(s,{\bf{R}})$ contains the rotation matrix ${\bf{R}}$. Therefore, the function $F(s,{\bf{R}})$ can be minimized by maximizing the item $s\sum\nolimits_{i=1}^{{N_{p}}}{\vec{d}_{i}^{T}}{\bf{R}}{\vec{m}_{i}}$. As $s$ is a scale factor, maximizing the item $s\sum\nolimits_{i=1}^{{N_{p}}}{\vec{d}_{i}^{T}}{\bf{R}}{\vec{m}_{i}}$ is equivalent to maximizing the item $\sum\nolimits_{i=1}^{{N_{p}}}{\vec{d}_{i}^{T}}{\bf{R}}{\vec{m}_{i}}$ for calculating the rotation matrix. According to [25], the matrix H with its SVD can be computed as follows:

[TABLE]

[TABLE]

where $svd(.)$ denotes the function of singular value decomposition. Then, the rotation matrix can be obtained as follows:

[TABLE]

(b) Scale factor

Taking the derivative of $F(s,{\bf{R}})$ with respective to the scale factor $s$, it is easy to get the following result:

[TABLE]

Let $\frac{{\partial F(s,{\bf{R}})}}{{\partial s}}=0$, the scale factor can be calculated as follows:

[TABLE]

Obviously, if the value of $n$ is set to be greater than 2 in Eq. (5), there are three items in ${\partial F(s,{\bf{R}})/\partial s}$. Accordingly, the scale factor $s$ will exit more than one solution, which may be the complex number. This is invalid and unexpected.

(c) Translation vector

After getting the rotation matrix and scale factor, the translation vector can be calculated by Eq. (8).

[TABLE]

By repeating Step(1)$\sim$(2), the scaling ICP algorithm can achieve the scaling registration of absolutely overlapping point sets.

3.2 The scaling and trimmed ICP algorithm

In many practical applications, point sets to be registered are always partially overlapping. Therefore, the proposed scaling registration algorithm should be extended to align partially overlapping point sets.

Suppose there are two overlapping overlapping point sets in $\mathbb{R}^{m}$, the model shape $Q\buildrel\Delta\over{=}\{{\vec{q}_{i}}\}_{i=1}^{{N_{q}}}$ and a data shape $P\buildrel\Delta\over{=}\{{\vec{p}_{i}}\}_{i=1}^{{N_{p}}}$ $({N_{p}},{N_{q}}\in\mathbb{N})$, where $\xi$ represents the overlapping percentage of the data shape, ${P_{\xi}}$ denotes the point subset, which is the overlapping part of the data shape to the model shape and $N_{p}^{{}^{\prime}}$ indicates the number of points in ${P_{\xi}}$. According to the idea of the trimmed ICP algorithm [9], the optimal scaling transformation can be obtained by minimizing the following objective function:

[TABLE]

where

[TABLE]

To solve Eq. (17), the scaling and trimmed ICP (sTrICP) algorithm can be proposed to achieve the scaling registration by iterations. Given the initial transformation $(s_{0},{\bf{R}_{0}},{\vec{t}_{0}})$, three steps are included in each iteration:

(1) Establish the point correspondence:

[TABLE]

(2) Update the overlapping percentage and its corresponding subset:

[TABLE]

(3) Calculate the scaling transformation:

[TABLE]

Obviously, Step (1) and (3) are similar to that of the proposed scaling registration algorithm. And Step (2) can be solved in a sequence processing manner. Denote ${\vec{p}_{i,k}}={s_{k-1}}{{\bf{R}}_{k-1}}{\vec{p}_{i}}+{\vec{t}_{k-1}}$, then the point pairs $\{{\vec{p}_{i,k}},{\vec{q}_{{c_{k}}(i)}}\}_{i=1}^{{N_{p}}}$ can be sorted by their distances in the ascending order. Each time, a pair of sorted points can be added to compute the corresponding value of $\Psi(\xi)$. By traveling all sorted point pairs, it is easy to obtain the minimum value $\Psi({\xi_{k}})$, which corresponds to the optimal overlapping percentage ${\xi_{k}}$. Then, the points involved in the front ${\xi_{k}}{N_{p}}$ sorted point pairs can be selected to update the corresponding point subset ${P_{{\xi_{k}}}}$.

By repeating Step(1)$\sim$(3), the scaling and trimmed ICP algorithm can achieve the scaling registration of partially overlapping point sets.

3.3 Proof of Convergence

As the sTrICP algorithm is proposed based on the idea of the original ICP algorithm, they have the similar convergence property. The following theorem will explain this convergence in detail so as to prove that this algorithm converges in theory.

Theorem. The sTrICP algorithm converges monotonically to a local minimum with respect to the objective function value.

The following proof shows that, for each step of the TsICP algorithm, the value of the objective function can be no worse than the previous step.

Proof. Suppose there are two partially overlapping point sets, $P\buildrel\Delta\over{=}\{{\vec{p}_{i}}\}_{i=1}^{{N_{p}}}$ and $Q\buildrel\Delta\over{=}\{{\vec{q}_{i}}\}_{i=1}^{{N_{q}}}$. Denote ${T_{k}}=({s_{k}},{{\bf{R}}_{k}},{\vec{t}_{k}})$, ${\xi_{k}}$ and ${P_{{\xi_{k}}}}$ as the scaling transformation, overlapping percentage, and overlapping part of $P$ to $Q$, respectively. In the first step of the $k$th iteration, the nearest neighbour ${\vec{q}_{{c_{k}}(i)}}$ in $Q$ is searched for the point ${\vec{p}_{i}}$ in $P$. Let ${\vec{p}_{i,k-1}}={s_{k-1}}{{\bf{R}}_{k-1}}{\vec{p}_{i}}+{\vec{t}_{k-1}}$. The value of objective function can then be defined as:

[TABLE]

In the second step of the $k$th iteration, the percentage ${\xi_{k}}$ and corresponding subset ${P_{{\xi_{k}}}}$ are updated. Then, the updated value of objective function is:

[TABLE]

As ${\xi_{k}}$ and ${P_{{\xi_{k}}}}$ are calculated from Eq. (20), it is reasonable to conclude that ${\eta_{k}}\leq{e_{k}}$. In the third step of the $k$th iteration, $\{{\vec{q}_{{c_{k}}(i)}}\}_{i=1}^{\left|{{P_{{\xi_{k}}}}}\right|}$ is registered with $\{{\vec{p}_{i,k-1}}\}_{i=1}^{\left|{{P_{{\xi_{k}}}}}\right|}$ and the scaling transformation ${{T}_{k}}{\rm{=(}}{s_{k}}{\rm{,}}{{\bf{R}}_{k}},{\vec{t}_{k}})$ is optimized. Let ${\vec{p}_{i,k}}={s_{k}}{{\bf{R}}_{k}}{\vec{p}_{i}}+{\vec{t}_{k}}$. And the value of objective function becomes:

[TABLE]

Since the scaling transformation ${{T}_{k}}{\rm{=(}}{s_{k}}{\rm{,}}{{\bf{R}}_{k}},{\vec{t}_{k}})$ is computed from Eq. (21), it is easy to conclude that ${\varepsilon_{k}}\leq{\eta_{k}}$.

In the $(k+1)$th iteration, the new nearest neighbour ${\vec{q}_{{c_{k+1}}(i)}}$ can be searched for ${\vec{p}_{i}}$ again. Subsequently, the new value of objective function can be defined as:

[TABLE]

Because $e_{k+1}$ is the result updated from the Eq. (19), therefore:

[TABLE]

Hence, repeating the above procedures, the following results can be reasonably obtained:

$0\leq...\leq{e_{k+1}}\leq{\varepsilon_{k}}\leq{\eta_{k}}\leq{e_{k}}\leq...$  for all $k$

According to the Monotonic Sequence Theorem: every bounded monotonic sequence of real numbers is convergent, we can reasonably conclude that the sTrICP algorithm always converges monotonically to a local minimum with respect to the objective function value.

4 Experimental results

In this section, experiments were performed on two public available datasets: (1) 2D shapes in part B of CE-Shape [27], (2) the Stanford 3D Scanning Repository [28].

To demonstrate its performance, the proposed approaches (sTrICP) was compared with the ICP algorithm with bounded scale [18] and the ICP algorithm based on the bidirectional distance measurement [20], which are abbreviated to TrICPBs and BiTrICP, respectively. The reason for the comparison with these two algorithms is that all these approaches are proposed under the framework of the original ICP algorithm. As the performance of TrICPBs is relatively affected by the boundaries of scale factor, the Principal Component Analysis (PCA) method was adopted to estimate its initial value $s_{0}$. For the comprehensive comparison, two groups of boundaries: $[0.9{s_{0}},1.1{s_{0}}]$ and $[0.5{s_{0}},2{s_{0}}]$ are both applied to the TrICPBs algorithm. For distinction, the TrICPBs algorithm with narrow and wide boundaries are abbreviated to TrICPBNs and TrICPBWs, respectively.

All the competing approaches utilized the nearest neighbor search method based on $k$-d tree to establish correspondences. Besides, no change of all point correspondences and the maximum iterations are both adopted as the termination condition of iterations for these completed approaches. Experiments were performed in MATLAB and conducted on a laptop with 2.5GHz processor of double-cores and 4GB RAM.

4.1 2D shapes

Here, experiments were conducted to illustrate the performance of the proposed approach for the scaling registration of 2D shapes. The data sets were selected from Part B of CE-Shape-1, which is a large 2D shapes database. Before the experiment, one shape was selected as the data shape and the other one was viewed as the model shape. Fig. 2 displays Apple shapes selected for the experiment. To generate partially overlapping point sets, some selected shapes should be cut by one part. Then, edge points were extracted for each shapes so as to obtain two point sets to be registered. Besides, a random scale factor was generated and might be applied to the data shape. Accordingly, all competing approaches can be applied to the registration of these shape pairs.

During experiment, the runtime required for the registration of each range image pair was recorded for all competing approaches. As these registration approaches are all efficient, 20 Monte Carlo (MC) trials were carried out for each shape pair so as to eliminate randomness. For comparison, Table 1 records the objective function value and average runtime for all these competing approaches.

As shown in Fig. 1, all competing approaches can obtain accurate results for the Apple shapes. Therefore, Fig. 3 illustrates the change of objective function value for all competing approaches tested on this shape. As depicted in Fig. 3, the sTrICP algorithm is locally convergent for the registration of 2D point sets. Besides, all these registration approaches converge in the similar way with a fast speed. This is because all these registration approaches are proposed under the framework of the original ICP algorithm, which can achieve registration with fast speed.

For other shapes, both the BiTrICP and sTrICP algorithms can also achieve the scaling registration with good accuracy. By introducing the bidirectional measurement, the BiTrICP algorithm can avoid the ill-posed problem involved in the scaling registration. Therefore, this approach is very robust. But it is somewhat time-consuming due to the establishment of the bidirectional correspondences. While, the sTrICP algorithm is proposed the new objective function, which imposes the emphasis on the scale factor. As shown in the Fig. 1, when the scale factor converges to 0, the mean square error gradually turns to be small and the value of the item $s^{2}$ reduces rapidly. Therefore, the objective function become large. Since the scaling transformation can be calculated by minimizing the well-designed objective function, the emphasis on the scale factor can avoid the ill-posed problem involved in the scaling registration. Obviously, there is no need to limit the boundaries of the scale factor. What’s more, it only required to establish the forward correspondence. Therefore, the proposed approach is also very robust and more efficient than that of the BiTrICP algorithm.

However, the TrICPBs algorithm may not always obtain good results for the scaling registration of 2D shape pairs. To avoid the ill-posed problem, the TrICPBs algorithm suggests that the boundaries of the scale factor should be limited due to its initial value. As displayed in the second row of Fig. LABEL:fig:2DRes, for partially overlapping point sets, it is difficult to estimate the initial value of the scale factor. To cover the real scale factor, the scale factor should be limited with a wide range. In this case, the scale factor also will converge to the direction of 0 and trapped into a local minimum around the low boundary. But if the scale factor is limited with a narrow range, which may not cover its real value, the TrICPBs approach will never obtain the accurate results for the scaling registration.

Therefore, the proposed approach can achieve the scaling registration of 2D shapes with superior performances than other two related approaches.

4.2 Range images

Here, experiments were conducted to illustrate the performance of the proposed approach for scaling registration of range images. The data sets were selected from the Stanford 3D Scanning Repository. Before experiment, one range image was selected as the data shape and the other one was viewed as the model shape. Besides, a random scale factor was generated and applied to the data shape. Subsequently, all competing approaches can be applied to registration of these range image pairs.

During experiment, the runtime required for the registration of each range image pair was recorded for all competing approaches. As all competing approaches are efficient, 10 MC trials were carried out for each range image pair so as to eliminate randomness. For comparison, Table 2 records the objective function value and average runtime for all these competing approaches.

As shown in Fig. 2, all competing approaches can obtain accurate results for the Stanford Happy Buddha. Therefore, Fig. 3 illustrates the change of objective function value for all all competing approaches tested on this data set. As depicted in Fig. 4, the sTrICP algorithm is locally convergent for the registration of 3D point sets. As mentioned before, the sTrICP algorithm shares the similar framework with other competing approaches. Therefore, they have almost the same convergent speed. Besides, both the sTrICP and BiTrICP algorithms can always obtain good scaling registration for other range image pairs. But the BiTrICP algorithm is required to establish bidirectional correspondences for each point in the data shape, it costs more time than the sTrICP algorithm, which only needs to establish the forward correspondence for each point in the data shape.

For the range image with high overlapping percentage, such as the Stanford Happy Buddha, the PCA method can provide good estimation of the scaler factor. Therefore, it is easy for the TrICPBs algorithm to obtain accurate results for scaling registration. For the range image with common overlapping percentages, such as the Stanford Bunny and Dragon, the PCA method is unable to provide good estimation of the scaler factor. Subsequently, a narrow range of boundaries may not cover the real scale factor, so the TrICPBs algorithm can never obtain the accurate results for scaling registration of range image pairs. Otherwise, a wide range of boundaries will make the TrICPBs algorithm easy to trap into local minimum around the low boundary. Therefore, the performance of the TrICPBs algorithm is more worse than that of other two approaches.

Therefore, the proposed approach can achieve the scaling registration of range images with the best performance among these competing approaches.

5 Application: Grid map merging

Mapping has been considered to be an important prerequisite of truly autonomous robots. Its efficiency can potentially be increased in a significant way by utilizing multiple robots exploring different parts of the environment. Therefore, a key question is how to merge maps built by the different robots into a single global map.

To verify the proposed approach, an experiment was performed on a real robot dataset, which is available for public access [30]. As shown in Fig. 5, there are two grid maps with different resolutions, which were produced by the robot exploring different parts of the same environment. To achieve grid maps, the scale-invariant feature transform (SIFT) features [13] were extracted and matched for two grid maps to be merged. Then, the random sample consensus (RANSAC) algorithm [26] was applied to confirm two geometrically consistent feature matches, which can be utilized to estimate the initial transformation for the scaling registration algorithm. As shown in Fig. 5, the initial transformation is inaccurate and the merging map is coarse. Subsequently, the edge points can be extracted from these two grid maps. According to our previous work [5], the problem of merging grid maps with different resolutions can be reasonably viewed as the scaling registration of the partially overlapping point sets. Therefore, the sTrICP algorithm can be applied to the registration of extracted edge point sets and obtain the accurate scaling transformation. Finally, these two grid maps can be merged into a single global map, which is more fine than the initial merged one.

To further illustrate the application, Fig. 6 also displays another grid map merging result, where the two grid maps to be merged are in different resolutions. After applying the sTrICP algorithm, we can obtain a global grid map, which has the similar resolution as the reference grid map.

6 Conclusions

This paper proposes a novel approach for the scaling registration of $m$-D point sets, which are absolutely or partially overlapping to each other. Firstly, for the absolutely overlapping point sets, a novel objective function is designed by imposing the emphasis on the scale factor. Similar to other registration problems, the new objective function also includes the item of the least-square error, which exits as its numerator. The significant difference is the introduction of the denominator item, which is the square of the scale factor and can avoid the ill-posed problem caused by the scale factor. Subsequently, a variant of the ICP algorithm is then proposed to solve this new objective function. What’s more, this approach is further extended to align partially overlapping point sets for the practical applications. Experimental results illustrate that the proposed approach can achieve the scaling registration of $m$-D point sets with good performance.

Acknowledgments

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