## 1.

## Introduction

Point set registration is important and complex in image processing. One robust and efficient approach to handle this problem is the iterative closest point (ICP) algorithm,^{1} but it cannot effectively register point sets with outliers including noise and missing data, which exist widely in applications.^{2} To overcome this problem, some scholars used thresholds or probability of distances to discard outliers,^{3} but these methods need an appropriate threshold which depends on the structure of point sets and thus is hard to choose. Meanwhile, some researchers attempted to reject outliers by adopting a coarse to fine process,^{4} but the methods often perform poorly in the case of a large amount of outliers in point sets.

Different from the above-mentioned work, Chetverikov
^{5} proposed the trimmed ICP (TrICP) algorithm for partially overlapping registration, which incorporates an overlapping percentage into a least square function to trim outliers. This algorithm, however, is time-consuming in that it needs to be repeated by traveling all possible overlapping percentage to search the best result even though Golden Section Search is used to accelerate the search speed. Based on this work, Phillips
^{6} presented the fractional ICP (FICP) algorithm for much faster speed, which computes the best correspondence and the overlapping percentage simultaneously, whereas the method depends greatly on a parameter which may lead to a false registration. To cope with the problem, this paper proposes a new objective function. Moreover, a novel robust and efficient ICP algorithm is presented to compute optimal overlapping percentage automatically and to complete the registration successfully.

## 2.

## Problem Statement

Given two point sets in
[TeX:]
\documentclass[12pt]{minimal}\begin{document}${\mathbb {R}^m}$\end{document}
${\mathbb{R}}^{m}$
, a model shape
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$M \buildrel\Delta \over = \lbrace \vec{m}_i \rbrace _{i = 1}^{N_m }$\end{document}
$M\genfrac{}{}{0pt}{}{\Delta}{=}{\left\{{\stackrel{\u20d7}{m}}_{i}\right\}}_{i=1}^{{N}_{m}}$
and a data shape
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$P \buildrel\Delta \over = \lbrace \vec{p}_i \rbrace _{i = 1}^{N_p } \quad (N_m ,N_p \in {\mathbb N})$\end{document}
$P\genfrac{}{}{0pt}{}{\Delta}{=}{\left\{{\stackrel{\u20d7}{p}}_{i}\right\}}_{i=1}^{{N}_{p}}\phantom{\rule{1em}{0ex}}({N}_{m},{N}_{p}\in \mathbb{N})$
. The rigid registration of *m*-D point sets is to build up correspondence and calculate rigid transformation (**R**, *t*) between *P* and *M*, hence

## 1

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \begin{array}{*{20}c} \mathop {\min }\limits _{{\bf {R}},\vec{t},c(i)\in \lbrace 1,2,\cdots ,N_m \rbrace } & {\displaystyle\sum \limits _{i = 1}^{N_p } {\Vert {({\bf {R}}\vec{p}_i + \vec{t}) - \vec{m}_{c(i)} } \Vert _2^2 } } \\[17pt] {s.t.} & {{\bf {R}}^{\rm {T}} {\bf {R}} = {\bf {I}}_m ,{\rm { }}\det ({\bf {R}}) = 1} \end{array} \end{equation} \end{document} $$\begin{array}{cc}\underset{\mathbf{R},\stackrel{\u20d7}{t},c\left(i\right)\in \{1,2,\cdots ,{N}_{m}\}}{min}& {\displaystyle \sum _{i=1}^{{N}_{p}}\Vert (\mathbf{R}{\stackrel{\u20d7}{p}}_{i}+\stackrel{\u20d7}{t})-{\stackrel{\u20d7}{m}}_{c\left(i\right)}{\Vert}_{2}^{2}}\\ s.t.& {\mathbf{R}}^{\mathrm{T}}\mathbf{R}={\mathbf{I}}_{m},\mathrm{det}\left(\mathbf{R}\right)=1\end{array}$$*c*(

*i*) is the index of the

*i*'th corresponding point in model shape

*M*.

For partially overlapping point sets *P* and *M*, assume *P*_{r} is the overlapping part of *P*, which can match with *M*, and *r* denotes the overlapping percentage of *P*. The goal of partial registration is to build up the correspondence of *P*_{r} and *M*, and the distance is
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\sum _{\vec{p}_i \in P_r } {\vert \vert {\rm {\bf R}}\vec{p}_i + \vec{t} - \vec{m}_{c(i)} \vert \vert _2^2 }$\end{document}
${\sum}_{{\stackrel{\u20d7}{p}}_{i}\in {P}_{r}}\left|\right|\mathbf{R}{\stackrel{\u20d7}{p}}_{i}+\stackrel{\u20d7}{t}-{\stackrel{\u20d7}{m}}_{c\left(i\right)}{\left|\right|}_{2}^{2}$
which increases with *r*. Hence, to include enough corresponding points as well as to exclude outliers, we divide the square distance by a penalty function *g*(*r*), an increasing function with respect to *r*. As too little overlapping part does not provide enough information, we assume the overlapping percentage *r* varies from 0.5 to 1. Therefore,

## 2

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \begin{array}{*{20}{l}} \mathop {\min }\limits _{c,r,{P_r},R,\vec{t} } \hfill & {\displaystyle\frac{1}{{{e^\lambda }{r^\lambda }}}\displaystyle\sum \limits _{{{\vec{p}}_i} \in {P_r}} {||{\mathbf {R}}{{\vec{p}}_i} + \vec{t} - {{\vec{m}}_{c(i)}}||_2^2} } \\[15pt] \; \; \; s.t. & {{{\mathbf {R}}^T}{\mathbf {R}} = {\mathbf {I}}_m,{\text{ }}\det ({\mathbf {R}}) = 1} \hfill \\[5pt] {} & {\lambda \in [{\lambda _{\min }},{\lambda _{\max }}],{\text{ }}r \in [0.5,1]} \hfill ,\\[5pt] {} & {{P_r} \in P,{\text{ }}|{P_r}| = r|P|} \hfill \\ \end{array} \end{equation} \end{document} $$\begin{array}{cc}\underset{c,r,{P}_{r},R,\stackrel{\u20d7}{t}}{min}\hfill & {\displaystyle \frac{1}{{e}^{\lambda}{r}^{\lambda}}\sum _{{\stackrel{\u20d7}{p}}_{i}\in {P}_{r}}\left|\right|\mathbf{R}{\stackrel{\u20d7}{p}}_{i}+\stackrel{\u20d7}{t}-{\stackrel{\u20d7}{m}}_{c\left(i\right)}{\left|\right|}_{2}^{2}}\\ \phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}s.t.\hfill & {\mathbf{R}}^{T}\mathbf{R}={\mathbf{I}}_{m},\phantom{\rule{4.pt}{0ex}}\mathrm{det}\left(\mathbf{R}\right)=1\hfill \\ & \lambda \in [{\lambda}_{\mathrm{min}},{\lambda}_{\mathrm{max}}],\phantom{\rule{4.pt}{0ex}}r\in [0.5,1],\\ & {P}_{r}\in P,\phantom{\rule{4.pt}{0ex}}|{P}_{r}|=r|P|\hfill \end{array}$$*e*

^{( · )}is the exponential function, λ is the control parameter, and | · | is the cardinality of a set, meaning the number of elements of the set.

In Eq. 2, it is easily proved that the objective function is a convex function with respect to *r*, which means the computation of *r* can be easily completed by iteration.

## 3.

## Point Set Registration with Outliers

ICP is good in minimizing Eq. 1 by repeating two steps: compute the closest point in model shape for each point in data shape and calculate the rigid transformation. As Eq. 2 is a convex function with respect to *r*, we compute *r* at each iterative step. Hence, similar to ICP, the objective function 2 can be optimized iteratively after the initialization:

**Step 1.** Set up the correspondence with the (*k* − 1)'th rigid transformation
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$({\rm {\bf R}}_{k - 1} ,\vec{t}_{k - 1} )$\end{document}
$({\mathbf{R}}_{k-1},{\stackrel{\u20d7}{t}}_{k-1})$
:

## 3

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} c_k (i) = \mathop {\arg \min }\limits _{j \in \lbrace 1,2, \ldots ,N_m \rbrace } \vert \vert {\rm {\bf R}}_{k - 1} \vec{p}_i + \vec{t}_{k - 1} - \vec{m}_j \vert \vert _2^2. \end{equation} \end{document} $${c}_{k}\left(i\right)=\underset{j\in \{1,2,...,{N}_{m}\}}{\mathrm{arg}\mathrm{min}}\left|\right|{\mathbf{R}}_{k-1}{\stackrel{\u20d7}{p}}_{i}+{\stackrel{\u20d7}{t}}_{k-1}-{\stackrel{\u20d7}{m}}_{j}{\left|\right|}_{2}^{2}.$$^{7}whose computation complexity is

*O*(

*N*

_{p}ln

*N*

_{m}).

**Step 2.** Compute the *k*'th overlapping percentage *r*_{k} and the subset
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$P_{r_k }$\end{document}
${P}_{{r}_{k}}$
according to the known control parameter λ and the corresponding point sets
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\lbrace {\rm {\bf R}}_{k - 1} \vec{p}_i + \vec{t}_{k - 1} \rbrace _{i = 1}^{N_p }$\end{document}
${\{{\mathbf{R}}_{k-1}{\stackrel{\u20d7}{p}}_{i}+{\stackrel{\u20d7}{t}}_{k-1}\}}_{i=1}^{{N}_{p}}$
and
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\lbrace \vec{m}_{c_k (i)} \rbrace _{i = 1}^{N_p }$\end{document}
${\left\{{\stackrel{\u20d7}{m}}_{{c}_{k}\left(i\right)}\right\}}_{i=1}^{{N}_{p}}$
by minimizing the objective functions:

## 4

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} r_k = \mathop {\arg \min }\limits _{0.5 \le r \le 1} \sum \limits _{\vec{p}_i \in P_r } {\vert \vert {\rm {\bf R}}_{k - 1} \vec{p}_i + \vec{t}_{k - 1} - \vec{m}_{c_k (i)} \vert \vert _2^2 } / (e^\lambda r^\lambda ), \end{equation} \end{document} $${r}_{k}=\underset{0.5\le r\le 1}{\mathrm{arg}\mathrm{min}}\sum _{{\stackrel{\u20d7}{p}}_{i}\in {P}_{r}}\left|\right|{\mathbf{R}}_{k-1}{\stackrel{\u20d7}{p}}_{i}+{\stackrel{\u20d7}{t}}_{k-1}-{\stackrel{\u20d7}{m}}_{{c}_{k}\left(i\right)}{\left|\right|}_{2}^{2}/\left({e}^{\lambda}{r}^{\lambda}\right),$$## 5

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} P_{r_k} = \mathop {\arg \min }\limits _{P_r \subset P,\vert P_r \vert = r \vert P\vert } \sum \limits _{\vec{p}_i \in P_r } {\vert \vert {\rm {\bf R}}_{k - 1} \vec{p}_i + \vec{t}_{k - 1} - \vec{m}_{c_k (i)} \vert \vert _2^2 }. \end{equation} \end{document} $${P}_{{r}_{k}}=\underset{{P}_{r}\subset P,|{P}_{r}|=r|P|}{\mathrm{arg}\mathrm{min}}\sum _{{\stackrel{\u20d7}{p}}_{i}\in {P}_{r}}\left|\right|{\mathbf{R}}_{k-1}{\stackrel{\u20d7}{p}}_{i}+{\stackrel{\u20d7}{t}}_{k-1}-{\stackrel{\u20d7}{m}}_{{c}_{k}\left(i\right)}{\left|\right|}_{2}^{2}.$$*P*

_{r}and compute the new corresponding

*f*(

*r*). Finally, we find the minimum

*f*(

*r*

_{k}) by traveling all

*f*(

*r*), and then select the first

*r*

_{k}

*N*

_{p}point pairs as the subset [TeX:] \documentclass[12pt]{minimal}\begin{document}$P_{r_k }$\end{document} ${P}_{{r}_{k}}$ .

**Step 3.** Compute the rigid transformation by minimizing the following squared distance:

## 6

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} ({{\bf {R}}^*},{\vec{t}^*}) {=}\!\! \mathop {\arg \min }\limits _{\scriptstyle {{\bf {R}}^{\rm {T}}}{\bf {R}} {=} {{\bf {I}}_m}, \scriptstyle \det ({\bf {R}}) {=} 1} \sum \limits _{{{\vec{p}}_i} \in {P_{{r_k}}}}\!\! {||{\bf {R}}({{\bf {R}}_{k {-} 1}}{{\vec{p}}_i} {+} {{\vec{t}}_{k - 1}}) {+} \vec{t} {-} {{\vec{m}}_{{c_k}(i)}}||_2^2},\!\!\!\!\! \nonumber\\ \end{eqnarray} \end{document} $$\begin{array}{c}\hfill ({\mathbf{R}}^{*},{\stackrel{\u20d7}{t}}^{*})=\underset{{\scriptstyle {\mathbf{R}}^{\mathrm{T}}\mathbf{R}={\mathbf{I}}_{m},\mathrm{det}\left(\mathbf{R}\right)=1}}{\mathrm{arg}\mathrm{min}}\sum _{{\stackrel{\u20d7}{p}}_{i}\in {P}_{{r}_{k}}}\left|\right|\mathbf{R}({\mathbf{R}}_{k-1}{\stackrel{\u20d7}{p}}_{i}+{\stackrel{\u20d7}{t}}_{k-1})+\stackrel{\u20d7}{t}-{\stackrel{\u20d7}{m}}_{{c}_{k}\left(i\right)}{\left|\right|}_{2}^{2},\end{array}$$## 7

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} {\rm {\bf R}}_k = {\rm {\bf R}}^\ast {\rm {\bf R}}_{k - 1} ,\mbox{ }\vec{t}_k = {\rm {\bf R}}^\ast \vec{t}_{k - 1} + \vec{t}^\ast . \end{equation} \end{document} $${\mathbf{R}}_{k}={\mathbf{R}}^{*}{\mathbf{R}}_{k-1},\phantom{\rule{4.pt}{0ex}}{\stackrel{\u20d7}{t}}_{k}={\mathbf{R}}^{*}{\stackrel{\u20d7}{t}}_{k-1}+{\stackrel{\u20d7}{t}}^{*}.$$^{8}

**Step 4.** Let
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\varepsilon _k = \sum _{\vec{p}_i \in P_{r_k } } {\vert \vert {\rm {\bf R}}_k \vec{p}_i + \vec{t}_k - \vec{m}_{c_k (i)} \vert \vert _2^2 } / (e^\lambda r_k^\lambda )$\end{document}
${\varepsilon}_{k}={\sum}_{{\stackrel{\u20d7}{p}}_{i}\in {P}_{{r}_{k}}}\left|\right|{\mathbf{R}}_{k}{\stackrel{\u20d7}{p}}_{i}+{\stackrel{\u20d7}{t}}_{k}-{\stackrel{\u20d7}{m}}_{{c}_{k}\left(i\right)}{\left|\right|}_{2}^{2}/\left({e}^{\lambda}{r}_{k}^{\lambda}\right)$
. If |ε_{k} − ε_{k − 1}| is sufficiently small or *k* reaches a maximum number of iterations, the inside loop stops with ϕ(λ) = ε_{k} and λ is updated by λ = λ − λ_{step}. Repeat steps 1 through 3 by decreasing the parameter λ from its maximal value λ_{max } to its minimal value λ_{min }.

Finally, ϕ(λ) decreases at the initial stage when λ increases from λ_{min } to λ_{max }, and the best registration result is located at the first increasing point of ϕ(λ).

**Discussion on the parameter λ.** For Eq. 2, *r* increases with λ, namely, a larger λ leads to more point pairs. Following this, we run the iteration by changing λ in a decreasing order, which keeps as many point pairs as possible at the beginning and filters outliers gradually. Otherwise, the algorithm may easily drop into a local minimum because of insufficient points.

For partial registration, using as many point pairs as possible will yield a more reliable and accurate result. When two point sets are matched, it is easily proved ϕ(λ) is a decreasing function with respect to λ. As λ varies from λ_{min } to λ_{max }, ϕ(λ) decreases at first and then increases when mismatching points are found. Therefore, the optimal result is located at the first increasing point of ϕ(λ) which guarantees enough information used in registration.

The proposed algorithm iterates in a similar way as the ICP algorithm. Despite more steps of our algorithm than the ICP algorithm at each iteration, both algorithms spend the same time on the correspondence computation which is the most time-consuming part. Compared with this, the calculation time of other steps can be relatively ignored. Therefore, our algorithm can get similar computation time as ICP, which is quite fast.

## 4.

## Experimental Results

We compare our algorithm with ICP,^{1} TrICP,^{5} and FICP,^{6} in which the latter two are given different values of λ, and the results are reported in root-mean-square (RMS) error.

## 4.1.

### 2D Shapes Registration

Experiments are conducted here on part B of CE-Shape-1 (Ref. 9), where random noise exists on shapes. First, the totally overlapping bat and horse shapes are used to demonstrate our algorithm is robust for full registration. Second, to show the robustness of our algorithm for partial registration, we select one shape pair and trim down one part of each shape to make them partially overlap. The compared results are displayed in Table 1.

## Table 1

Compared results on two-dimensional shapes.

bat | horse | trimmed bat | trimmed horse | |||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

Algorithms | λ | r | RMS | λ | r | RMS | λ | r | RMS | λ | r | RMS |

ICP | – | 1 | 0.55 | – | 1 | 0.52 | – | 1 | 19.41 | – | 1 | 13.69 |

TrICP | 2 | 0.98 | 0.51 | 2 | 0.51 | 6.94 | 2 | 0.69 | 0.49 | 2 | 0.88 | 6.44 |

TrICP | 4 | 0.98 | 0.51 | 4 | 0.64 | 9.99 | 4 | 0.69 | 0.49 | 4 | 0.88 | 6.44 |

TrICP | 6 | 0.98 | 0.51 | 6 | 0.71 | 13.52 | 6 | 0.69 | 0.49 | 6 | 0.88 | 6.44 |

FICP | 1.3 | 0.5 | 9.11 | 1.3 | 0.5 | 6.80 | 1.3 | 0.51 | 15.02 | 1.3 | 0.5 | 8.39 |

FICP | 3 | 1 | 0.55 | 3 | 0.70 | 12.59 | 3 | 0.70 | 0.51 | 3 | 1 | 14.26 |

FICP | 5 | 1 | 0.55 | 5 | 1 | 0.52 | 5 | 1 | 19.41 | 5 | 1 | 13.69 |

Ours | 6 | 1 | 0.53 | 6 | 1 | 0.51 | 6 | 0.70 | 0.50 | 4 | 0.80 | 0.59 |

As is shown in Table 1, our method always yields a good registration with small RMS. ICP only deals with registration of totally overlapping shapes and performs poorly for trimmed shapes that have a large amount of missing data. TrICP can complete some partial registration though, it is quite unstable for different shapes. As TrICP searches the best overlapping percentage by traveling all possible values, the results are likely to be similar. FICP achieves on some shapes, but it substantially relies on the preset value of λ which varies for different shapes registration. Therefore, both TrICP and FICP cannot complete full and partial registration well. In addition, our method automatically obtains an appropriate value of parameter λ, which not only guarantees a proper overlapping percentage *r* but also uses as many point pairs as possible to obtain more accurate registration results. In bat and horse shapes registration, the largest overlapping percentage reaches while RMS is small. One intutive example of the registration results is shown in Fig. 1.

## 4.2.

### 3D Range Data Registration

We compare our algorithm with ICP and FICP on Stanford 3D Scanning Repository^{10} which is acquired from different viewpoints in practice. As the database has too many points and TrICP is too slow and unstable, we do not give the results of TrICP. In the database, Stanford Bunny (bun000 with 40,256 points and bun 045 with 40,097 points), Dragon (dragonStandRight_0 with 41,841 points and dragonStandRight_24 with 34,836 points), and Happy Budda (happyStandRight_0 with 78,056 points and happyStandRight_24 with 75,582 points) are adopted. The compared results are shown in Table 2.

## Table 2

Compared results on Stanford Database.

Bunny | Dragon | Happy Budda | |||||||
---|---|---|---|---|---|---|---|---|---|

Algorithms | λ | r | RMS (× 10−3) | λ | r | RMS (× 10−3) | λ | r | RMS (× 10−3) |

ICP | – | 1 | 2.05 | – | 1 | 1.86 | – | 1 | 2.28 |

FICP | 0.95 | 0.50 | 0.87 | 0.95 | 0.50 | 0.60 | 0.95 | 0.50 | 0.69 |

FICP | 3 | 0.91 | 0.38 | 3 | 0.90 | 0.35 | 3 | 0.82 | 0.26 |

FICP | 5 | 0.93 | 0.40 | 5 | 0.92 | 0.38 | 5 | 0.83 | 0.29 |

Ours | 5 | 0.91 | 0.35 | 5 | 0.90 | 0.32 | 5 | 0.81 | 0.24 |

Table 2 shows that our method gives the best registration results with least RMS, and computes the parameter λ automatically while ICP produces worse results and FICP performs well only when λ is equal to 3 or 5, which demonstrates FICP depends on the parameter λ and its results are unstable. The registration results are displayed in Fig. 2.

## 5.

## Conclusion

This paper proposes a novel approach for registration between two *m*-D point sets with outliers in the way of incorporating an overlapping percentage into the ICP algorithm. Compared with previous works, our algorithm can automatically compute rigid transformation, correspondence, and overlapping percentage without influence of a parameter. Furthermore, the best overlapping percentage is computed by including as much information as possible, which leads to a satisfying registration. A series of compared experiments demonstrate that our algorithm is stable and precise.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant Nos. 61005014, 90920008, and 90920301.

## References

## Biography

**Shaoyi Du** received double Bachelor degrees in computational mathematics and in computer science in 2002 and received his MS degree in applied mathematics in 2005 and PhD degree in pattern recognition and intelligence system from Xi'an Jiaotong University, China in 2009. His research interests include computer vision, pattern recognition and machine learning.

**Jihua Zhu** received his BE degree from the Central South University China in 2004 and PhD degree in control science and engineering from Xi'an Jiaotong University, China, in 2011. He is a lecturer of Xi'an Jiaotong University. His research interests include mobile robot and image registration.

**Nanning Zheng** graduated in 1975 from the Department of Electrical Engineering, Xi'an Jiaotong University, China, and received his ME degree in information and control engineering from Xi'an Jiaotong University, China in 1981, and PhD degree in electrical engineering from Keio University, Japan, in 1985. He is currently a professor and the director of the Institute of Artificial Intelligence and Robotics in Xi'an Jiaotong University. His research interests include computer vision, pattern recognition, computational intelligence, image processing, and hardware implementation of intelligent systems.

**Yuehu Liu** received his BE and ME degrees from Xi'an Jiaotong University, China in 1984 and 1989 respectively, and PhD degree in electrical engineering from Keio University, Japan, in 2000. He is a professor of Xi'an Jiaotong University. His research interests include computer vision, computer graphics and digital content analysis.

**Ce Li** graduated in 1995 from School of Electrical Engineering and Automation of Hefei University of Technology, China, and received his MS degree in information and control engineering from Xi'an Jiaotong University, China, in 2003. He is currently a PhD student in the institute of Artificial Intelligence and Robotics in Xi'an Jiaotong University. His research interests include computer vision, pattern recognition.