# 3DMNDT: 3D multi-view registration method based on the normal distributions
transform

Jihua Zhu*1, Di Wang1, Jiaxi Mu2, Huimin Lu3, Zhiqiang Tian1, Zhongyu Li1 *
Corresponding author: Jihua Zhu (email: zhujh@xjtu.edu.cn). 1 School of
Software Engineering, Xi’an Jiaotong University, China
2 School of Artificial Intelligence, Xi’an Jiaotong University, China
2 Kyushu Institute of Technology, Japan

###### Abstract

The normal distributions transform (NDT) is an effective paradigm for the
point set registration. This method is originally designed for pair-wise
registration and it will suffer from great challenges when applied to multi-
view registration. Under the NDT framework, this paper proposes a novel multi-
view registration method, named 3D multi-view registration based on the normal
distributions transform (3DMNDT), which integrates the K-means clustering and
Lie algebra solver to achieve multi-view registration. More specifically, the
multi-view registration is first cast into the problem of maximum likelihood
estimation. Then, the K-means algorithm is utilized to divide all data points
into different clusters, where a normal distribution is computed to locally
models the probability of measuring a data point in each cluster.
Subsequently, the registration problem is formulated by the NDT-based
likelihood function. To maximize this likelihood function, the Lie algebra
solver is developed to sequentially optimize each rigid transformation. The
proposed method alternately implements data point clustering, NDT computing,
and likelihood maximization until desired registration results are obtained.
Experimental results tested on benchmark data sets illustrate that the
proposed method can achieve state-of-the-art performance for multi-view
registration.

###### Index Terms:

Multi-view registration, normal distributions transform, Lie algebra, K-means
clustering.

##  I Introduction

Point set registration is a fundamental problem arising from main domains,
such as computer vision [1, 2, 3], robotics [4, 5, 6], and medical image
analysis [7, 8]. With the development of scanning technology, many scan
devices are able to acquire accurate 3D scan data from a scene or object.
However, they are difficult to scan the whole object or scene at one time due
to the limited field of view. For 3D reconstruction, it is necessary to
acquire scan data at different viewpoints and align these scans, which arises
the problem of point set registration. Given the reference frame, the goal of
multi-view registration is to compute the spatial transformations of each
point set, so as to unify them into the same coordinate system.

In the literature, Iterative Closest Point (ICP) [9] is one of the most
popular solutions to point set registration. Given initial rigid
transformation, this method alternatively establishes point correspondence and
optimizes rigid transformation to achieve pair-wise registration with good
accuracy. However, the original ICP method is unable to register partially
overlapping point sets and it is time-consuming due to the registration over
raw data points. Then, the normal distributions transform (NDT) [10, 11] is
proposed to solve pair-wise registration problem. For the pair-wise
registration, it divides the space occupied by the target point set into
regularized cells and approximates data points in each cell by the normal
distribution. Further, each data point of source point set is assigned to one
NDT and then the registration problem is formulated by the NDT-based objective
function, which is optimized by the Newton’s solver. Since one NDT is utilized
to approximate all data points in each cell, this method is efficient to align
large scale data sets. What’s more, it is able to align partially overlapping
point sets. Therefore, it is very popular in computer vision.

Although NDT has been widely applied to pair-wise registration, it has not
been extended to 3D multi-view registration, which suffers from great
challenges. Firstly, it is difficult to directly formulate the multi-view
registration problem and simultaneously optimize a lot of rigid
transformations under the original NDT framework. Secondly, the original NDT
divides the space into regularized cells. This regularized division leads to
the large variance of the number of points in each cell, which reduces the
registration performance of NDT. Thirdly, the original NDT utilizes the
Newton’s solver, which is nasty to directly compute Jacobian matrix and
Hessian matrix for the optimization of 3D rigid transformations.

In this paper, we extend the original NDT and propose a novel method named
3DMNDT to solve the 3D multi-view registration problem. The contributions of
this paper are delivered as follows: 1) We formulate the multi-view
registration problem into NDT based likelihood function. Therefore, multi-view
registration can be achieved by the maximum likelihood estimation. 2) We
utilize the K-means clustering to replace the regularized division of space
cubes, which can avoid the unbalanced division of data points. 3) We propose
the Lie algebra solver to sequentially optimize each rigid transformation.
Compared with the Newton’s solver, this optimization algorithm avoids the
calculation of Jacobian matrix and Hessian matrix. 4) We test 3DMNDT on seven
benchmark data sets and compare it with several state-of-the-art methods.
Experimental results illustrate its superior performance for 3D multi-view
registration.

The remainder of this paper is organized as follows. Section 2 briefly
discusses related methods of point set registration. In Section 3, we
formulate the multi-view registration problem under the NDT framework and then
derive the 3DMNDT to solve 3D multi-view registration. We also analyze its
computation complexity and present its implementation details. Section 5
displays the experimental results on benchmark data sets. Finally, we conclude
this paper in Section 6.

##  II Related work

This section only discusses these methods related to the proposed method for
multi-view registration. For convenience, the terms motion and rigid
transformation are utilized interchangeably throughout this paper.

According to the number of point sets being registered, point set registration
can be divided into pair-wise registration problem and multi-view registration
problem. For pair-wise registration, many variants have been derived from the
original ICP algorithm and they improve the performance of ICP from different
aspects [12]. For partially overlapping point sets, the Trimmed ICP (TrICP)
algorithm [13] introduces the overlap parameter to automatically distinguish
overlapping regions and is able to achieve good registration. For accurate
registration, Myronenko and Song proposed the Coherent Point Drift (CPD)
algorithm [14], which considers pair-wise registration as a probability
density estimation problem. More specifically, it takes the source point set
as the Gaussian mixture model (GMM) centroids and fits them to the target
point set by maximizing the likelihood. Besides, Jian et al. [15] formulated
the pair-wise registration problem as aligning two GMMs, where each point set
is represented by one GMM. Then it achieves pair-wise registration by
minimizing $L_{2}$ measure between two GMMs. Similar to CPD, Gao and Tedrake
[16] proposed the FilterReg, which takes the target point set as GMM centroids
and fits the source point set to GMM by maximizing the likelihood. Although
these GMM-based methods are powerful for both rigid and non-rigid
registration, it requires to build all point correspondences between two point
sets being registered, which leads to high computational complexity.

Although these methods are effective for pair-wise registration, they are
local convergent and can only obtain desired registration results in the case
of provided with good initial parameters. To estimate initial parameters,
Sipiran and Bustos [17] proposed the 3D Harris detection method, which can
detect the interest-point from 3D point set represented by meshes. Given
detected interest points, feature descriptors with match method [18, 19, 20]
can be applied to achieve pair-wise registration for point sets with arbitrary
orientations. Recently, some deep learning methods [21, 22, 23] have been
proposed to achieve global pair-wise registration. Although some of these
methods are effective, most of them are supervised learning method, which
requires many data sets to train the registration model. What’s more, these
methods cannot be directly extended to solve multi-view registration problem.

Usually, multi-view registration is much more complex than pair-wise
registration but receives less attention. For multi-view registration, the
primary method [24] is to incrementally register and merge point set pairs
until all point sets are merged into one model. Although this method is easy
to operate, it suffers from the accumulation error problem with the increase
of merged point sets. Then, Bergevin et al. [25] proposed the first method for
multi-view registration. This method sequentially optimizes the rigid
transformation over the established point correspondences, which are built
between one point set between each other point set. This method overcomes the
problem of cumulative error, but it neglects non-overlapping regions between
each point set pairs, which reduces the registration accuracy. To address this
issue, Zhu et al. [26] proposed the coarse-to-fine registration method, which
sequentially optimizes the trimmed errors of point correspondences, which
built between one point set and other aligned point sets. With the increase of
point sets, the coarse-to-fine registration method tends to be trapped into
local minimum. Accordingly, Tang et al. [27] proposed hierarchical multi-view
registration method, which continuously implements multi-view registration on
small number of point sets and merges them as the one in the next
implementation.

As there are many rigid transformations involved in multi-view registration,
most of methods utilize alternating optimization strategy to sequentially
optimize each rigid transformation. However, some methods utilize the batch
optimization strategy. Given point correspondences established between each
point set pair, Krishnan et al. [28] proposed a multi-view registration method
via optimization-on-a-manifold, which can simultaneously optimize all rigid
transformations. However, it is difficult to obtain point correspondences
between each point set pair, which should be provided in advance for this
method. Accordingly, Mateo et al. [29] treated pairwise correspondences as
missing data and formulated multi-view registration as the maximum likelihood
estimation, where all rigid transformations are optimized by the
Expectation–Maximization (EM) algorithm in batch mode. Since there are massive
point correspondences, this method requires to estimate massive hidden
variables, which is time-consuming.

Since multi-view registration is more difficult than pair-wise registration,
one feasible method is to recover the global motions from relative motions,
where the term motion denotes rigid transformation. Accordingly, Govindu et
al. [30] proposed the motion averaging algorithm and applied it to multi-view
registration. This method utilizes the Lie group structure of motions [31] to
implement the averaging of all available relative motions. Theoretically,
accurate global motions can be obtained from the number of relative motions,
which is equal to or more than the number of point sets being registered.
However, this method is sensitive to outliers, where even one outlier will
lead to the failure of multi-view registration. Further, Arrigoni et al. [32,
33] cast multi-view registration to a low-rank and sparse (LRS) matrix
decomposition problem, where registration results are recovered from the low-
rank matrix. In this method, it requires to concatenate all available relative
motions as well as some zero matrices into one large matrix, which is then
decomposed into a low-rank matrix and a sparse matrix. This method is robust
to outliers but it tends to be failed when the concatenated matrix is sparse
[33]. Besides, these two methods equally treat all relative motions, which may
possess different reliabilities. Accordingly, Guo et al. [34] proposed weighed
motion averaging algorithm and Jin et al. [35] proposed the weighted LRS for
multi-view registration. Since these two methods pay more attention to these
relative motions with high reliability, they are more likely to obtain desired
registration results than their original methods.

Recently, Georgios et al. [36] cast multi-view registration to the clustering
problem and proposed JRMPC to solve this special problem. This method assumes
that all data points are generated from a central GMM and utilizes the EM
algorithm to estimate the GMM components as well as rigid transformations.
Since this method requires to estimate massive model parameters related to all
data points, it is time-consuming and likely to be trapped into local minimum.
To address these issues, Zhu et al. [37] extended the k-means algorithm to
solve multi-view registration problem. Compared with JRMPC, this method is
efficient and more likely to obtain desired registration results. However, it
utilizes one cluster centroid to approximate all data points in the same
cluster, which inevitably leads to a lot of information loss and reduces the
registration accuracy. Further, Zhu et al. [38] proposed an efficient GMM-
based method under the perspective of EM. This method assumes that each data
point is drawn from its corresponding GMM, where its nearest neighbors ( NNs )
in other point sets are taken as Gaussian centroids with similar covariance
and membership probability. Since this method only requires to estimate quite
a few rigid transformations as well as one covariance, it is efficient than
other mixture model-based methods. However, these mixture model-based methods
are more time-consuming than most previous methods.

The original NDT was first proposed to solve the registration of 2D point sets
in robotics [10]. As this method utilizes the ND to approximate a set of point
sets, it is efficient to achieve pair-wise registration. Further, it was
extended to align 3D point sets [11], where the computation of Jacobian and
Hessian matrix with respect to Euler angles is really nasty. Besides, the
spatial division mechanism leads to regular cells, which reduces the
registration performance. Accordingly, Arun and Steven [39] replaced regular
division by k-means clustering to divide all data points of target point set
into different clusters, where the target point set is always static during
registration. However, this method cannot be directly extended to multi-view
registration, where rigid transformations of point set are constantly changed
and all data points should be dynamically clustered at each iteration. In this
paper, we formulate multi-view registration problem into the NDT-based
objective function, which is then maximized by the integration of the k-means
clustering and the Lie algebra solver.

##  III Problem formulation

Let ${\bf
V}_{i}=[{{\bf{v}}_{i,1}}...{{\bf{v}}_{i,j}}...{{\bf{v}}_{i,{N_{i}}}}]$
indicates $N_{i}$ data points in the $i$-th point set and
${\mathbb{V}}=\\{{{\bf V}_{i}}\\}_{i=1}^{M}$ denotes $M$ point sets being
registered. Similar to some previous methods, we regard multi-view
registration as the clustering problem. Given the reference frame, the goal of
multi-view registration is to simultaneously divide all data points into $K$
clusters and estimate accurate rigid transformations ${\mathbb{T}}=\\{{\bf
R}_{i}\in{\rm SO(3)},{\bf t}_{i}\in{\mathbb{R}}^{3}\\}_{i=1}^{M}$ for each
point set. For convenience, the rigid transformation is concatenated into
motion as follows:

${\cal T}_{i}=\begin{bmatrix}{\bf{R}}_{i}&{\bf{t}}_{i}\\\ {\bf
0}&1\end{bmatrix},$ (1)

where ${\cal T}_{i}\in{\rm SE(3)}$. Besides, we define the group action of
${\rm SE(3)}$ on a 3D vector ${\bf v}_{i,j}$ as ${\cal T}_{i}{\bf
v}_{i,j}={\bf R}_{i}{\bf v}_{i,j}+{\bf t}_{i}$. For each cluster, it is
convenient to compute the NDT $({\bm{\mu}}_{k},{\bm{\Omega}}_{k})$, which
denote the mean (centroid) and information matrix (the inverse of covariance
matrix), respectively.

Suppose the aligned data point ${{\cal T}_{i}}{{\bf{v}}_{i,j}}$ belongs the
$c_{(i,j)}$-th cluster, then it is reasonable to assume that the aligned data
point ${{\cal T}_{i}}{{\bf{v}}_{i,j}}$ is drawn from the normal distribution
${\cal T}_{i}{\bf v}_{i,j}\sim{\cal
N}({\bm{\mu}}_{c(i,j)},{\bm{\Omega}}_{c(i,j)}^{-1})$, where $c_{(i,j)}$
denotes the index of cluster, which contains the data point
${{\bf{v}}_{i,j}}$. Therefore, the probability function is defined for
${{\bf{v}}_{i,j}}$ as follows:

$p({\bf v}_{i,j})={\cal N}({\cal T}_{i}{\bf
v}_{i,j};{\bm{\mu}}_{c(i,j)},{\bm{\Omega}}^{-1}_{c(i,j)}).$ (2)

Further, all data points are assumed to be independent of each other and their
joint probability distribution is formulated as:

$p({\mathbb{V}})=\prod_{i=1}^{M}\prod_{j=1}^{N_{i}}{\cal N}({\cal T}_{i}{\bf
v}_{i,j};{\bm{\mu}}_{c(i,j)},{\bm{\Omega}}_{c(i,j)}^{-1}),$ (3)

where the normal distribution on ${\mathbb{R}}^{3}$ is defined as:

${\cal N}({\bf
x};{\bm{\mu}},{\bm{\Omega}}^{-1})=\sqrt{\frac{|{\bm{\Omega}}|}{(2\pi)^{3}}}\exp\\{-\frac{1}{2}({\bf
x}-{\bm{\mu}})^{\top}{\bm{\Omega}}({\bf x}-{\bm{\mu}})\\}.$

This joint probability distribution is utilized to define the logarithm of
likelihood function as:

$\displaystyle L({\bm{\Theta}})$ $\displaystyle=\log p({\mathbb{V}})$ (4)
$\displaystyle=\sum_{i=1}^{M}\sum_{j=1}^{N_{i}}\log{\cal N}({\cal T}_{i}{\bf
v}_{c(i,j)};{\bm{\mu}}_{c(i,j)},{\bm{\Omega}}_{c(i,j)}^{-1}),$

where ${\bm{\Theta}}\buildrel\Delta\over{=}\left\\{{\left({\cal
T}_{i}\right)_{i=1}^{M},}{\left({{{\bm{\mu}}_{k}},{{\bm{\Omega}}_{k}}}\right)_{k=1}^{K}}\right\\}$
denotes the model parameters.

Accordingly, the multi-view registration is transformed into the maximum
likelihood estimation, which is formulated as:

$\displaystyle\max_{\bm{\Theta}}L({\bm{\Theta}})$ (5) $\displaystyle{\rm
s.t.}\quad{\cal T}_{i}\in{\rm SE(3)},i=1,2,\cdots,M.$

By the optimization of Eq. (5), multi-view registration can be achieved with
data point clustering.

Figure 1: Illustration of the proposed method, which integrates the K-means
clustering and Lie algebra solver to achieve multi-view registration under the
NDT framework. Given initial model parameters, the proposed method uses the
K-means clustering algorithm to divide all aligned data points into different
clusters and compute the NDT for each cluster. Then it utilizes the Lie
algebra solver to sequentially optimize each rigid transformation. The desired
registration results will be obtained by alternatively implementing data point
clustering, NDT computation, and rigid transformation optimization.

##  IV 3DMNDT method

For the optimization of Eq. (5), we propose the 3D multi-view version of NDT
method (3DMNDT) displayed in Fig. 1. As shown in Fig. 1, this method optimizes
Eq. (5) by iterations. Given initial rigid transformations $({\cal
T}^{0}_{i})_{i=1}^{M}$ and initial cluster centroids
$({\bm{\mu}}_{k}^{0})_{k=1}^{K}$ , three steps are alternately executed in
each iteration.

(1) Divide all data points into different clusters.

$c^{h}(i,j)=\mathop{\arg\min}_{k=1,2,\cdots,K}{\left\|{\cal T}_{i}^{h-1}{\bf
v}_{i,j}-{\bm{\mu}}_{k}^{h-1}\right\|_{2}^{2}}.$ (6)

(2) Compute the NDT for each cluster.

${\bm{\mu}}_{k}^{h}=\frac{\sum_{i=1}^{M}\sum_{j=1}^{N_{i}}\\{c^{h}(i,j)=k\\}{\cal
T}_{i}^{h-1}{\bf
v}_{i,j}}{\sum_{i=1}^{M}\sum_{j=1}^{N_{i}}\\{c^{h}(i,j)=k\\}},$ (7)
${\bm{\Omega}}_{k}^{h}=({\bm{\Sigma}}_{k}^{h}+\varepsilon{\bf I}_{3})^{-1},$
(8)

where $\varepsilon=10^{-6}$ prevents the occurrence of zero denominator and
${\bm{\Sigma}}^{h}_{k}$ indicates the covariance matrix calculated by:

${\bm{\Sigma}}_{k}^{h}=\frac{\sum_{i=1}^{M}\sum_{j=1}^{N_{i}}\\{c^{h}(i,j)=k\\}\Delta{\bf
v}_{i,j}^{h}(\Delta{\bf
v}_{i,j}^{h})^{\top}}{\sum_{i=1}^{M}\sum_{j=1}^{N_{i}}\\{c^{h}(i,j)=k\\}},$
(9)

where $\Delta{\bf v}_{i,j}^{h}={\cal T}_{i}^{h-1}{\bf
v}_{i,j}-{\bm{\mu}}_{k}^{h}$.

It should be noted that some clusters may only contain a small number of data
points. In 3D cases, points of one cluster may be perfectly co-planar, the
covariance matrix will always be singular and can not be inverted. For this
reason, NDTs only computed for clusters containing more than 5 points.
Meanwhile, other clusters are regarded as invalid clusters.

(3) Rigid transformation optimization.

If one data point belongs to a invalid cluster, it should be neglected in the
optimization of rigid transformation. Accordingly, all rigid transformations
can be optimized by maximum likelihood estimation:

$\max_{\mathbb{T}}\sum\limits_{i=1}^{M}\sum_{j=1}^{N_{i}^{{}^{\prime}}}\log{\cal
N}({\cal T}_{i}{\bf
v}_{i,j};{\bm{\mu}}_{c^{h}(i,j)},{\bm{\Omega}}^{-1}_{c^{h}(i,j)}),$ (10)

where $N^{{}^{\prime}}_{i}\leq N_{i}$ denotes the number of valid data point
in the $i$-th point set.

To obtain the desired results, the proposed method should iteratively operate
these three steps until the likelihood $L({\bm{\Theta}})$ has no significant
change or the iteration number $h$ reaches the maximum value. Obviously, Steps
1 denotes the nearest neighbor search problem, which can be efficiently solved
by the $k$-d tree based method [40]. Step 2 is a trivial calculation problem.
However, Step 3 remains to be a difficult optimization problem, which is the
most critical step of the proposed method.

###  IV-A Optimizing by the Lie algebra solver

It seems difficult to simultaneously optimize all rigid transformations
involved in Eq. (5). Given clustering results, each rigid transformation can
be separately optimized by the following objective function:

$\max_{{\cal T}_{i}\in{\rm SE(3)}}\sum_{j=1}^{N_{i}^{{}^{\prime}}}\log{\cal
N}({\cal T}_{i}{\bf
v}_{i,j};{\bm{\mu}}_{c^{h}(i,j)},{\bm{\Omega}}^{-1}_{c^{h}(i,j)}),$ (11)

where $i=1,2,\cdots,M$. For simplicity, we discard the iteration number $h$,
the point set number $i$ and cluster number $c(i,j)$ in all variables.
Accordingly, the above optimization can be simplified as a non-linear least-
squares problem:

$\displaystyle\max_{{\cal T}\in{\rm SE(3)}}L({\cal T})$
$\displaystyle=\max_{{\cal T}\in{\rm SE(3)}}\sum_{j=1}^{N}\log{\cal N}({\cal
T}{\bf v}_{j};{\bm{\mu}}_{j},{\bm{\Omega}}^{-1}_{j})$ (12)
$\displaystyle=\min_{{\cal T}\in{\rm SE(3)}}\sum_{j=1}^{N}{\bf
r}_{j}^{\top}{\bm{\Omega}}_{j}{\bf r}_{j},$

where ${\bf r}_{j}={\cal T}{\bf v}_{j}-{\bm{\mu}}_{j}$ indicates the residual
error between the data point ${\cal T}{\bf v}_{j}$ and its cluster centroid
${\bm{\mu}}_{j}$.

To circumvent the constraints stems from rigid transformation ${\cal T}$, one
feasible method is to use Lie algebra as follows:

$\max_{{\cal T}\in{\rm SE(3)}}L({\cal
T})=\max_{{\bm{\xi}}\in{\mathbb{R}}^{6}}L(\exp({\bm{\xi}}^{\wedge}){\cal
T}^{0}),$ (13)

where ${\cal T}^{0}$ denotes rigid transformation estimated from previous
iteration, ${\bm{\xi}}$ denotes the vector space of Lie algebra
${\mathfrak{s}e(3)}$,
$(\cdot)^{\wedge}:{\mathbb{R}}^{6}\rightarrow{\mathbb{R}}^{4\times 4}$ denotes
the skew operator which converts vector space ${\mathbb{R}}^{6}$ to Lie
algebra ${\mathfrak{s}e(3)}$, $\exp(\cdot):{\mathfrak{s}e(3)}\rightarrow{\rm
SE(3)}$ denotes the exponential map which lifts Lie algebra to Lie Group. As
the exponential map of ${\mathfrak{s}e(3)}$ involves the calculation of
Jacobian matrix, it is quite complex and highly non-linear. Accordingly, we
turn to utilize retraction technique [41].

In the retraction technique, the vector space of ${\mathfrak{s}e(3)}$ is
decomposed into rotational and translational components, i.e.,
${\bm{\xi}}=\\{{\bm{\xi}}_{R}\in{\mathbb{R}}^{3},{\bm{\xi}}_{t}\in{\mathbb{R}}^{3}\\}$,
and then the retraction technique ${\cal
R}(\cdot):{\mathbb{R}}^{6}\rightarrow{\rm SE(3)}$ lifts ${\bm{\xi}}$ to ${\rm
SE(3)}$ as follows:

${\cal
R}({\bm{\xi}})=\begin{bmatrix}\exp({\bm{\xi}}_{R}^{\wedge})&{\bm{\xi}}_{t}\\\
{\bf 0}&1\end{bmatrix},$ (14)

where the exponential map $\exp(\cdot):{\rm so(3)}\rightarrow{\rm SO(3)}$
lifts ${\mathfrak{s}o}(3)$ to ${\rm SO(3)}$. It is worth noting that the skew
operator has an important property:

${\bm{\xi}}_{R}^{\wedge}{\bf v}=-{\bf v}^{\wedge}{\bm{\xi}}_{R},\forall{\bf
v}\in{\mathbb{R}}^{3},$ (15)

which will be useful later on.

Since the group action of ${\rm SE(3)}$ is matrix multiplication, we denote
the current rigid transformation as ${\cal T}^{0}$, and assume the
perturbation of rigid transformation is ${\bm{\xi}}\in{\mathbb{R}}^{6}$.
Consequently, the optimization is with respect to ${\bm{\xi}}$, and the rigid
transformation is updated as follows:

${\cal T}^{*}={\cal R}({\bm{\xi}}^{*}){\cal T}^{0},$ (16)

where ${\bm{\xi}}^{*}$ denotes the optimized perturbation.

Subsequently, we present the following theoretical result. Proposition 1 By
utilizing retraction technique in Eq. (14), the optimal perturbation can be
attained by solving the following quadratic program (QP) problem:

$\displaystyle{\bm{\xi}}^{*}$
$\displaystyle=\mathop{\arg\min}\limits_{{\bm{\xi}}\in{\mathbb{R}}^{6}}{\bm{\xi}}^{\top}{\bf
H}{\bm{\xi}}+2{\bf b}^{\top}{\bm{\xi}}+c$ (17) $\displaystyle{\bf H}$
$\displaystyle=\sum_{j=1}^{N}{\bf H}_{j}^{\top}{\bm{\Omega}}_{j}{\bf
H}_{j}\in{\mathbb{R}}^{6\times 6}$ $\displaystyle{\bf b}$
$\displaystyle=\sum_{j=1}^{N}{\bf H}_{j}^{\top}{\bm{\Omega}}_{j}{\bf
r}_{j}^{0}\in{\mathbb{R}}^{6}$ $\displaystyle c$
$\displaystyle=\sum_{j=1}^{N}({\bf r}_{j}^{0})^{\top}{\bm{\Omega}}_{j}{\bf
r}_{j}^{0}\in{\mathbb{R}}^{+}$

where ${\bf H}_{j}\buildrel\Delta\over{=}\left[-\left({\bf R}^{0}{\bf
v}_{j}+{\bf t}^{0}\right)^{\wedge}\quad{\bf I}\right]\in{\mathbb{R}}^{3\times
6}$ and ${\bf r}_{j}^{0}={\cal T}^{0}{\bf v}_{j}-{\bm{\mu}}_{j}$.

Proof 1 Based on the assumption of small perturbation, the exponential
operator of ${\mathfrak{s}o}(3)$ can be approximated as:

$\exp\left({\bm{\xi}}_{R}^{\wedge}\right)\approx{\bf
I}+{\bm{\xi}}_{R}^{\wedge}.$ (18)

We can substitute the Eq. (18) into the Eq. (16), which leads to the rigid
transformation defined as:

$\displaystyle{\bf R}$ $\displaystyle=\exp({\bm{\xi}}_{R}^{\wedge}){\bf
R}^{0}\approx{\bf R}^{0}+{\bm{\xi}}_{R}^{\wedge}{\bf R}^{0}$ (19)
$\displaystyle{\bf t}$ $\displaystyle=\exp({\bm{\xi}}_{R}^{\wedge}){\bf
t}^{0}+{\bm{\xi}}_{t}\approx({\bf I}+{\bm{\xi}}_{R}^{\wedge}){\bf
t}^{0}+{\bm{\xi}}_{t}.$

Then, the residual error is formulated as:

$\displaystyle{\bf r}_{j}$ $\displaystyle={\cal T}{\bf v}_{j}-{\bm{\mu}}_{j}$
(20) $\displaystyle\approx({\bf R}^{0}+{\bm{\xi}}_{R}^{\wedge}{\bf R}^{0}){\bf
v}_{j}+{\bf t}^{0}+{\bm{\xi}}_{R}^{\wedge}{\bf
t}^{0}+{\bm{\xi}}_{t}-{\bm{\mu}}_{j}$
$\displaystyle={\bm{\xi}}_{R}^{\wedge}({\bf R}^{0}{\bf
v}_{j})+{\bm{\xi}}_{R}^{\wedge}{\bf t}^{0}+{\bm{\xi}}_{t}+{\bf r}_{j}^{0}$

Given the property depicted in Eq. (15), the residual error is converted into:

$\displaystyle{{\bf{r}}_{j}}=\begin{bmatrix}-({\bf R}^{0}{\bf v}_{j}+{\bf
t}^{0})^{\wedge}&{\bf I}\end{bmatrix}{\bm{\xi}}+{\bf r}_{j}^{0}={\bf
H}_{i}{\bm{\xi}}+{\bf r}_{j}^{0}.$ (21)

By replacing ${{\bf{r}}_{j}}$ with term ${\bf H}_{i}{\bm{\xi}}+{\bf
r}_{j}^{0}$, Eq. (12) can be converted into a quadratic program (QP) problem
depicted in Eq. (17), which is easy to be solved.

Obviously, the objective function of QP problem is a local approximation to
the original objective function, thus ${\bf H}$ is Hessian matrix and
${\bf{b}}$ is the gradient vector. The convergence of above optimization
behaves like a second-order optimization, which is faster than typical
gradient-based methods. And the optimal perturbation is calculated as follows:

${\bm{\xi}}^{*}=-{\bf H}^{\dagger}{\bf b},$ (22)

where $(\cdot)^{\dagger}$ denotes the pseudo-inverse operator. Finally, the
rigid transformation can be updated by Eq. (16), which is the optimal solution
of Eq. (12). For each point set, the Lie algebra solver can be sequentially
utilized to optimize its rigid transformation.

###  IV-B Implementation

As a local convergent method, 3DMNDT requires to be provided with initial
rigid transformations. Given data points with arbitrary orientation, we
utilize the coarse registration method [42] to estimate initial rigid
transformation for multi-view registration. Different from some previous
methods, this method should be additionally initialized with $K$ cluster
centroids. Therefore, we uniformly sample $K$ data points from all initial
aligned point sets and take them as the initial centroids to start iteration.
To obtain desired results, the proposed method should alternately operate
cluster update, ND computation, and transformation optimization until all
rigid transformations have no significant change or the iteration number $h$
reaches the maximum value $H$.

Based on the above description, the proposed method is summarized in Algorithm
1, where we set $H=300$. In the proposed method, there is a free parameter
$K$, which requires to be predefined. Frankly speaking, it is difficult to set
its optimal value. However, we can set the averaging number of data points
contained in each cluster, which can indirectly determine the cluster number.
Specifically, we set $(6+M)$ to be the averaging number of data points due to
two reasons. First, this setting can ensure that most clusters contain at
least 6 data points and reduce the invalid clusters, even when the number of
point sets is small. Second, the cluster number remains stable and will not
increase with the increase of the point set number.

Algorithm 1 3DMNDT method

0: Point sets ${\mathbb{V}}=\\{{{\bf V}_{i}}\\}_{i=1}^{M}$, initial
transformations ${\mathbb{T}}^{0}$.

1: Uniformly initialize ${\bm{\mu}}^{0}_{k}$, set $h=0,K=M+6$;

2: repeat

3: $h=h+1$;

4: Update cluster $c^{h}(i,j)$ by Eq. (6);

5: Compute NDT $({\bm{\mu}}_{k}^{h},{\bm{\Omega}}_{k}^{h})$ by Eqs. (7) and
(8).

6: for $\left(i=2:M\right)$ do

7: Optimize the perturbation ${{\bm{\xi}}_{i}^{*,h}}$ by Eq. (22);

8: Update the transformation ${\cal T}_{i}^{h}$ by Eq (16);

9: end for

10: Compute likelihood $L({\bm{\Theta}}^{h})$ by Eq. (4).

11: until
($\left|{L({{\bm{\Theta}}^{h}}){\rm{-}}L({{\bm{\Theta}}^{h-1}})}\right|<\varepsilon$)
or $\left(h>H\right)$

TABLE I: The total computation complexity of the proposed method in each operation Operation | Complexity
---|---
Build $k$-d tree | $O(HK\lg K)$
Build correspondence | $O(HN\lg K)$
Compute NDT | $O(HN)$
Optimize transformation | $O(HN)$

###  IV-C Computation complexity

This section analyzes the computation complexity of 3DMNDT. Before analysis,
we restate that there are $M$ point sets being registered, where the $i$-th
point set contains $N_{i}$ data points and the total data points is
$N=\Sigma_{i=1}^{M}N_{i}$. Besides, these data points desire to be divided
into $K$ clusters and the number of valid clusters is $K^{{}^{\prime}}<K$. For
the convenience of analysis, we assume the iteration number of this method is
$H$. In each iteration, three operations are included to optimize one rigid
transformation.

Update cluster. Before updating cluster, it requires to build $k$-d tree for
$K$ cluster centroids, which leads to the total complexity of $O(HK\lg K)$ for
$H$ iterations. Then, each data point should be assigned to one cluster by the
NN search, which leads to the total complexity of $O(HN\lg K)$ for $H$
iterations.

Compute NDT. After updating cluster, the NDT should be computed for each
cluster. Since each data point is only assigned to one cluster and each
cluster utilizes all its data points to compute one NDT, the complexity is
$O(N)$ in each iteration. For $H$ iterations, the total complexity is $O(HN)$

Optimize rigid transformation. The proposed method utilizes
$N_{i}^{{}^{\prime}}\approx N_{i}$ point-centroid pairs to optimize the $i$th
rigid transformation. To estimate $M$ rigid transformations, the total
complexity is $O(HN)$ for $H$ iterations.

Accordingly, Table I displays the total computation complexity of the proposed
method. As shown in Table I, the highest complexity operation is linear
proportion to $N$, $H$ and $lgK$. Usually, the cluster number $K$ is much
smaller than the point number $N$. What’s more, $\lg K$ is much smaller than
the cluster number $K$ . Therefore, the proposed method is efficient.

##  V Experiment results

In this section, we evaluate the proposed 3DMNDT method on various benchmark
data sets, which include: (1) Object data sets. Four data sets are downloaded
from the Stanford 3D Scanning Repository and two data sets are provided by
Torsello. Each of them was acquired from one object model in multiple view-
points. (2) Environment data set. Gazebo data set (Courtesy of Francois
Pomerleau) is acquired from outdoor environment and it is provided for robot
simultaneous localization and mapping (SLAM). Table II displays the scan
number and unit for each data set. In addition to scans, the ground truth of
rigid transformations are also provided to evaluate the performance of
registration method. For efficiency, each data set is uniformly down-sampled
to around 2000 points per scan.

To illustrate its performance, we compare 3DMNDT with a number of baselines:

(1) MATrICP [30]. This method utilizes the motion averaging algorithm to
recover global motions from a set of relative motions, which are estimated by
the TrICP algorithm.

(2) JRMPC [36]. This method supposes that all data points are drawn from a
central GMM and utilizes the EM algorithm to estimate the GMM components as
well as rigid transformations.

(3) TMM [36]. This method supposes that all data points are drawn from a
central Student’s t-mixture model (TMM) and utilizes the EM algorithm to
estimate the TMM components as well as rigid transformations.

(4) LRS [33] This method concatenates available relative motions into one
matrix, which is decomposed by the low rank and sparse (LRS) decomposition
algorithm to recover global motions for multi-view registration.

(5) K-means [37]. This method casts multi-view registration into clustering
problem and extends the k-means clustering algorithm to achieve multi-view
registration.

(6) EMPMR [38]. This method assumes that each data point is generated by one
unique GMM, where all Gaussian centroids are efficiently searched from each
other scans by the $k$-d based method.

TABLE II: Details of bench mark data sets utilized in the experiment. Dataset | Angel | Armadillo | Bunny | Buddha | Dragon | Hand | Gazebo
---|---|---|---|---|---|---|---
Scans | 36 | 12 | 10 | 15 | 15 | 36 | 32
Unit | dm | mm | mm | mm | mm | dm | m

Experimental results are reported in the form of rotation error and
translation error defined as:

$e_{\bf R}=\frac{1}{M}\sum_{i=1}^{M}{\cos^{-1}({\rm Tr}({\bf R}_{m,i}{\bf
R}_{g,i}^{\top})-1)\over 2},$ (23) $e_{t}=\frac{1}{M}\sum_{i=1}^{M}\|{\bf
t}_{m,i}-{\bf t}_{g,i}\|_{2},$ (24)

where ${\rm Tr}(\cdot)$ denotes the trace operator on a square matrix, ${\rm
cos}^{-1}:[-1,1]\rightarrow[0,\pi]$ denotes the inverse cosine function,
$\\{{\bf R}_{g,i},{\bf t}_{g,i}\\}$ and $\\{{\bf R}_{m,i},{\bf t}_{m,i}\\}$
denote the ground truth and the estimated one of the $i$-th rigid
transformation, respectively. Therefore, the unit of $e_{\bf R}$ is radian and
the unit of $e_{t}$ is similar to the unit of each data set. All compared
methods are implemented on Matlab without any extra library and utilize the
$k$-d tree method to search the NN. Experiments are carried out on a quad-core
3.6 GHz computer with 8 GB of memory.

TABLE III: Error comparison of different NDT variants, where bold numbers denote the best performance of each data set. Method |  | Angel | Armadillo | Bunny | Buddha | Dragon | Hand | Gazebo
---|---|---|---|---|---|---|---|---
Initial | $e_{\bf R}$ | 0.0221 | 0.0234 | 0.0239 | 0.0262 | 0.0251 | 0.0282 | 0.0241
$e_{t}$ | 2.0388 | 2.5333 | 2.1260 | 1.6535 | 1.5216 | 0.4945 | 0.1069
NDTO | $e_{\bf R}$ | 0.0258 | 0.0057 | 0.0113 | 0.0232 | 0.0349 | 0.0214 | 0.0422
$e_{t}$ | 1.4885 | 2.3317 | 1.0811 | 4.7053 | 6.5298 | 0.6159 | 0.5948
NDTL | $e_{\bf R}$ | 0.0422 | 0.0054 | 0.0041 | 0.0095 | 0.0084 | 0.0174 | 0.0650
$e_{t}$ | 2.1972 | 0.8756 | 0.6865 | 0.9500 | 1.6755 | 1.7976 | 0.3092
3DMNDT | $e_{\bf R}$ | 0.0008 | 0.0049 | 0.0024 | 0.0073 | 0.0103 | 0.0014 | 0.0082
$e_{t}$ | 0.1287 | 0.6688 | 0.2426 | 1.0712 | 1.3531 | 0.0876 | 0.0188

###  V-A Validation of clustering

To validate the proposed 3DMNDT, we compared it with two variants of NDT
algorithm, where NDTO sequentially registers each scan to other aligned scans
by original NDT and NDTL replaces the Newton-based solver with the Lie algebra
based solver in NDTO to achieve registration. Experimental results are
reported in the form of registration errors, which are displayed in Table III.
As shown in Table III, 3DMNDT can obtain more accurate registration results
among three competed methods.

Figure 2: Illustration of NDTs for different methods, where cyan points denote
these data points belong to invalid clusters. (a) Results of NDTO. (b) Result
of NDTL. (c) Result of 3DMNDT. Figure 3:  Convergence results of 3DNDT tested
on object data sets. For view convenience, we display the log-likelihood value
subtracted by the initial value for each data set.

Usually, data points are non-uniformly distributed in the occupied space. As
shown in Fig. 2, different division strategies lead to different NDT results.
In NDTO and NDTL, the occupied space are divided into regular cells, which
inevitably lead to unbalanced division. Specifically, some cells may only
contain a small number of data points, i.e. smaller than 6. To avoid singular
covariance, data points located in these cells should be neglected to compute
NDT and optimize rigid transformations, which inevitably reduces the
performance of multi-view registration. Accordingly, both NDTO and NDTL are
difficult to achieve promising results for multi-view registration. What’s
more, the Newton-based solver in NDTO is difficult to be implemented and
debugged due to two reasons: Firstly, the Jacobian and Hessian matrices
frequently generate many pages of dense algebra that should be converted into
code, which is a tedious, difficult and error-prone process. Secondly, the
Newton-based solver also involves complex line search, which requires hundreds
lines of code.

While, the k-means clustering is able to avoid unbalanced division and
guarantee that most clusters contain a certain number of points. Since 3DMNDT
utilizes k-means clustering, most of data points are utilized to compute NDTs
and optimize rigid transformations. Compared with the Newton’s solver, the Lie
algebra solver is a more direct and explicit mechanism for the optimization of
3D rigid transformation. Fig. 3 also displays the convergence of 3DMNDT with
respect to iterations. As shown in Fig. 3, the log-likelihood value increases
with the increase of iterations. Thanks for K-means clustering and the Lie
algebra solver, which enable 3DMNDT to achieve promising registration results
within acceptable iterations.

###  V-B Object data sets

To illustrate its registration performance, 3DMNDT is tested on object data
sets and compared with a number of baselines on accuracy, efficiency, and
robustness.

TABLE IV: Error comparison of all competed methods, where bold numbers denote the best performance of each data set. Method |  | Angel | Armadillo | Bunny | Buddha | Dragon | Hand
---|---|---|---|---|---|---|---
Initial | $e_{\bf R}$ | 0.0221 | 0.0234 | 0.0239 | 0.0262 | 0.0251 | 0.0282
$e_{t}$ | 2.0388 | 2.5333 | 2.1260 | 1.6535 | 1.5216 | 0.4945
MATrICP | $e_{\bf R}$ | 0.0103 | 0.0237 | 0.0088 | 0.0085 | 0.0116 | 0.0262
$e_{t}$ | 2.8651 | 4.2915 | 0.7181 | 0.8581 | 1.1037 | 1.3390
JRMPC | $e_{\bf R}$ | 0.0065 | 0.0121 | 0.0166 | 0.0161 | 0.0154 | 0.0044
$e_{t}$ | 1.9816 | 3.0651 | 1.7571 | 0.9629 | 1.7466 | 0.7069
TMM | $e_{\bf R}$ | 0.0062 | 0.0061 | 0.0080 | 0.0106 | 0.0154 | 0.0072
$e_{t}$ | 2.5197 | 1.9646 | 0.7031 | 1.2137 | 1.0689 | 0.7682
LRS | $e_{\bf R}$ | 0.0087 | 0.0410 | 0.0125 | 0.0181 | 0.0162 | 0.0265
$e_{t}$ | 2.0633 | 9.4194 | 0.9586 | 1.0597 | 1.5358 | 1.3825
K-means | $e_{\bf R}$ | 0.0056 | 0.0066 | 0.0100 | 0.0147 | 0.0134 | 0.0059
$e_{t}$ | 0.9789 | 2.3207 | 1.6088 | 1.0276 | 1.6140 | 0.4615
EMPMR | $e_{\bf R}$ | 0.0012 | 0.0052 | 0.0035 | 0.0067 | 0.0090 | 0.0014
$e_{t}$ | 0.2161 | 0.7145 | 0.3439 | 1.1023 | 1.1433 | 0.0887
3DMNDT | $e_{\bf R}$ | 0.0008 | 0.0049 | 0.0024 | 0.0073 | 0.0103 | 0.0015
$e_{t}$ | 0.1287 | 0.6688 | 0.2426 | 1.0712 | 1.3531 | 0.0876
Figure 4: Multi-view registration results in the form of cross-section, where
corresponding regions are illustrated on the aligned 3D models and thin curves
denote good registration results. (a) Aligned 3D models. (b) Initial results.
(c) Results MATrICP. (d) Results of JRMPC. (e) Results of TMM. (f) Results of
LRS. (g) Results of K-means. (h) Results of EMPMR. (i) Results of 3DMNDT.

####  V-B1 Accuracy

To compare registration accuracy, these registration methods are tested on six
object data sets, experimental results are reported in registration errors.
Table IV illustrates registration errors of all compared methods tested on six
object data sets. For more intuitive comparison, Fig. 4 displays all
registration results in the form of cross-section. As shown in Table IV and
Fig. 4, 3DMNDT and EMPMR are able to obtain the most promising registration
results for most of data sets.

For the multi-view registration, both MATrICP and LRS recover global motions
from a set of available relative motions, which are estimated by the TrICP
algorithm. Given accurate relative motions, these two methods may obtain
promising registration results. Since the TrICP algorithm is not error free,
relative motions inevitably contain errors, which reduce the accuracy of
multi-view registration. What’s more, MATrICP is sensitive to unreliable
relative motions, even one unreliable relative motions will lead to its
registration failure, e.g. the Armadillo data set. Meanwhile, LRS is dependent
on the ratio of available relative motions and low ratio will lead to
registration failure, e.g. the Armadillo data set. As probabilistic methods,
both JRMPC and TMM are expected to obtain very accurate registration results.
However, this is not the case. Because these two methods require to estimate
massive mixture model parameters as well as rigid transformations, which makes
them to be easily trapped into local minimum. While, EMPMR only requires to
estimate rigid transformations as well as one Gaussian variance, so it is more
likely to obtain promising registration results.

Although both K-means and 3DMNDT utilize k-means clustering to achieve multi-
view registration, their performance is totally different. In K-means, all
data points in each cluster are approximated by one centroid, which inevitably
leads to information loss. As K-means directly aligns each point set to all
cluster centroids for registration, it is difficult to obtain promising
results. While, 3DMNDT computes NDT for each cluster and achieves multi-view
registration by maximizing NDT-based likelihood function. Since NDT contains
statistical information of each cluster, 3DMNDT is more likely to obtain
accurate registration results.

####  V-B2 Efficiency

To compare registration efficiency, these registration methods are tested on
six object data sets. Experimental results are reported in the form of
averaging run time over 30 independent tests. Fig. 5 illustrates the
efficiency comparison of these registration methods tested on object data
sets. As shown in Fig. 5, K-means is the most efficient one among all
registration methods and 3DMNDT is comparable with K-means especially for
large data sets. Besides, TMM and JRPMC are more time-consuming than other
registration methods.

Figure 5: Efficiency comparison over 10 independent tests. For view
convenience, we display relative run time for all competed methods, where TMM
is most time-consuming method with (17.2717, 5.4772, 6.5912, 9.276, 10.4413,
15.9024) minutes corresponding to 100% run time of each data set.

Figure 6: Comparison of averaging rotation errors for each object data set.
(a) Angle. (b) Armadillo. (c) Bunny. (d) Buddha. (e) Dragon. (f) Hand.

Figure 7: Comparison of averaging translation errors for each object data set.
(a) Angle. (b) Armadillo. (c) Bunny. (d) Buddha. (e) Dragon. (f) Hand.

Both TMM and JRPMC assume that all data points are generated from a central
mixture model. To achieve multi-view registration, EM algorithm is utilized to
estimate all mixture components as well as rigid transformations. Since each
mixture component and rigid transformation are related to all data points, the
estimation is very time-consuming. Compared with JRPMC, TMM requires to
estimate another set of hidden variables. Therefore, TMM is the most time-
consuming one among all compared methods. As aforementioned, both MATrICP and
LRS recover global motions from available relative motions. Therefore, it
requires to build point correspondences between most scan pairs. Since the
convergence speed of MATrICP is slower than that of LRS, it is less efficient
than LRS. With the increase of scan number, MATrICP even becomes less
efficient than JRPMC. Besides, MATrICP is less efficient than EMPMR. Although
EMPMR requires to build point correspondences between each scan pair, its
convergence speed is faster than that of MATrICP. Accordingly, it is more
efficient than MATrICP but a little less efficient than LRS.

Similar to K-means, 3DMNDT also utilizes K-means clustering algorithm to
achieve multi-view registration. Specifically, K-means compute a centroid and
minimize the sum function of point-centroid errors to optimize rigid
transformations. While, 3DMNDT requires to compute NDT and minimize the NDT-
based likelihood function. Therefore, 3DMNDT is a little less efficient than
K-means, but is more efficient than other registration methods, especially for
large data sets.

####  V-B3 Robust to initialization

To illustrate its robustness to initialization, 3DMNDT as well as other
registration methods are also tested on object data sets with different
initial rigid transformations, where rotation angles are added with three
levels of uniformly distributed noises. Considering randomness, each group of
experiment is carried out by 20 independent tests. Experimental results are
reported in the form of averaging registration errors. Figs. 6 and 7
demonstrate averaging rotation errors and translation error of all competed
methods, respectively. Considering both rotation and translation errors,
3DMNDT is generally the most robust one to initialization among all competed
methods.

Both MATrICP and LRS recover global motions from a set of relative motions,
which are estimated by the pair-wise registration algorithm. Given a set of
reliable relative motions, MATrICP can obtain promising registration results.
However, this method is sensitive to unreliable relative motions. Without good
initialization, the probability of obtaining unreliable relative motions will
increase, which will lead to the failure of MATrICP. Although LRS is robust to
unreliable relative motions, it is sensitive to the ratio of reliable relative
motions. Under poor initialization, the number of reliable relative motions
will be reduced, which will lead to the failure of LRS. Therefore, LRS is also
sensitive to initialization. Since JRMPC and TMM requires to estimate massive
model parameters, they are more likely trapped into local minimum, especially
when the initialization is poor. Therefore, both of them are also sensitive to
initialization.

Benefiting from data clustering, both K-means and 3DMNDT have wide convergent
domains. However, K-means is difficult to obtain promising registration
results due to information loss. Since 3DMNDT utilizes one NDT to represent
data points of each cluster, it is more likely to obtain promising results due
to less information loss. Thereby, 3DMNDT is the most robust one to
initialization among all competed methods.

TABLE V: Comparison of registration error (mean) for competed methods tested on six data sets perturbed by two level of Gaussian noises, where bold numbers denote the best performance.  Method | Error | 50dB | 25dB
---|---|---|---
Angel | Armadillo | Bunny | Buddha | Dragon | Hand | Angel | Armadillo | Bunny | Buddha | Dragon | Hand
MATrICP | $e_{\bf R}$ | 0.0102 | 0.0204 | 0.0089 | 0.0086 | 0.0114 | 0.0261 | 0.0103 | 0.0204 | 0.0119 | 0.0123 | 0.0162 | 0.0369
$e_{t}$ | 2.8710 | 4.0146 | 0.6996 | 0.8424 | 1.0464 | 1.3417 | 2.8726 | 3.8764 | 0.7047 | 0.8433 | 1.0534 | 1.3419
JRMPC | $e_{\bf R}$ | 0.0062 | 0.0119 | 0.0146 | 0.0162 | 0.0153 | 0.0044 | 0.0062 | 0.0118 | 0.0165 | 0.0154 | 0.0157 | 0.0048
$e_{t}$ | 1.9069 | 3.0249 | 1.7402 | 0.9690 | 1.7823 | 0.7075 | 1.9628 | 3.0543 | 1.7404 | 0.9694 | 1.7923 | 0.7178
TMM | $e_{\bf R}$ | 0.0061 | 0.0060 | 0.0084 | 0.0102 | 0.0153 | 0.0072 | 0.0062 | 0.0062 | 0.0086 | 0.0104 | 0.0155 | 0.0073
$e_{t}$ | 2.4665 | 1.9844 | 0.7487 | 1.1904 | 1.0662 | 0.7792 | 2.4865 | 1.9913 | 0.7626 | 1.1950 | 1.0686 | 0.7802
LRS | $e_{\bf R}$ | 0.0086 | 0.0397 | 0.0127 | 0.0180 | 0.0150 | 0.0265 | 0.0123 | 0.0398 | 0.0127 | 0.0181 | 0.0153 | 0.0265
$e_{t}$ | 2.0599 | 9.1893 | 0.9882 | 1.0443 | 1.4367 | 1.3803 | 2.0559 | 9.1627 | 0.9905 | 1.0551 | 1.5105 | 1.3817
K-means | $e_{\bf R}$ | 0.0056 | 0.0064 | 0.0104 | 0.0140 | 0.0141 | 0.0055 | 0.0057 | 0.0067 | 0.0107 | 0.0142 | 0.0151 | 0.0056
$e_{t}$ | 1.0117 | 2.3448 | 1.6936 | 1.0684 | 1.3913 | 0.4674 | 1.0750 | 2.4698 | 1.7043 | 1.1212 | 1.5158 | 0.4690
EMPMR | $e_{\bf R}$ | 0.0012 | 0.0052 | 0.0035 | 0.0068 | 0.0090 | 0.0014 | 0.0012 | 0.0055 | 0.0037 | 0.0069 | 0.0090 | 0.0020
$e_{t}$ | 0.2191 | 0.7170 | 0.3439 | 1.0972 | 1.1179 | 0.0839 | 0.2253 | 0.7202 | 0.3444 | 1.1001 | 1.1443 | 0.1426
3DMNDT | $e_{\bf R}$ | 0.0007 | 0.0046 | 0.0025 | 0.0076 | 0.0108 | 0.0014 | 0.0008 | 0.0046 | 0.0027 | 0.0079 | 0.0109 | 0.0026
$e_{t}$ | 0.1258 | 0.6551 | 0.2449 | 1.0899 | 1.3507 | 0.0838 | 0.1488 | 0.6643 | 0.2519 | 1.0956 | 1.3729 | 0.1567
TABLE VI: Comparison of localization accuracy for different methods tested on Gazebo data set, where bold numbers denote the best performance. | Initial | K-means | JRMPC | LRS | MATrICP | TMM | EMPRM | 3DMNDT
---|---|---|---|---|---|---|---|---
$e_{\bf R}$ (rad.) | 0.0241 | 0.0097 | 0.0216 | 0.0092 | 0.0089 | 0.0145 | 0.0087 | 0.0082
$e_{t}$ (m) | 0.1069 | 0.0318 | 0.1005 | 0.0557 | 0.0511 | 0.0630 | 0.0185 | 0.0188
Time (min.) | / | 0.8817 | 7.8113 | 1.6324 | 1.9533 | 19.2284 | 0.4258 | 0.6832

####  V-B4 Robustness to data noise

To illustrate its robustness to data noises, 3DMNDT as well as other
registration methods are tested on object data sets added with random Gaussian
noises. Considering randomness, each group of experiment is carried out by 30
independent tests. Experimental results are reported in the form of averaging
registration errors. Tables V demonstrate comparison results under two levels
of noises. As illustrated in Tables V, the accuracy of all registration
methods decreases with the increase of noise level. In general, EMPMR and
3DMNDT are more robust to data noise than other registration methods.

Theoretically, JRMPC and TMM are probabilistic methods, where data noises are
considered in the mixture model. However, they are more likely trapped into
local minimum due to massive model parameters required to be optimized.
Therefore, they seem to be sensitive to data noise. Compared with JRMPC and
TMM, EMPMR requires to optimize fewer model parameters and it is more likely
to obtain promising registration results. Since EMPRM directly formulates data
noise in the objective function, it is very robust to data noise. Although
both MATrICP and LRS do not take data noise into consideration, they seem
robust to data noise for some data sets. This is because these two methods
recover global motions from a set of relative motions and the influence of
data noise on multi-view registration may be eliminated by the motion
averaging or LRS matrix decomposition algorithm. However, high level noises
may result in unreliable relative motions, which can lead to the failure of
multi-view registration. Meanwhile, K-means also does not consider data
noises, so it is less robust than EMPMR.

Different from EMPRM, 3DMNDT does not directly formulate data noises in the
objective function. However, it utilizes NDT, which allows data to contain
noise. Therefore, it is also robust to data noise and is able to achieve
promising registration results under different levels of Gaussian noises.

###  V-C Environment data set

Usually, multi-view registration methods can be potentially applied to SLAM.
To illustrate its SLAM performance, 3DMNDT as well as other methods are tested
on Gazebo data set. This data set is acquired by mobile robot equipped with
the laser range finder in outdoor environment, where the robot mainly moved on
the 2D plane ground and tracked the path to form a closed loop ($4\times
5\times 0.09$m). During robot movement, 31 range scans were acquired with
ground truth of location, which were recorded by IMU and GPS.

Figure 8: Illustration of localization results for different methods tested on
Gazebo data set. (a) Initialization. (b) K-means result. (c) LRS result. (d)
MATrICP result. (e) TMM result. (f) JRMPC result. (g) EMPRM result. (h) 3DMNDT
result.

Before registration, data point indicated ground surface should be filtered by
the height value and initial rigid transformations are provided by the TrICP
algorithm. Table VI illustrates the localization errors of all compared
methods. For visual comparison, Fig. 8 displays corresponding localization
results for these methods. As shown in Table VI and Fig. 8, 3DMNDT as well as
EMPRM is more accurate and efficient than other compared methods. It is worth
noting that both MATrICP and LRS only tend to reduce the localization error
around the closed loop. This is because the closed loop introduces more
available relative motions, which tend to substantially reduce the
localization error around the closed loop. While, other methods are able to
reduce localization error along the whole robot path. Among these methods,
3DMNDT and EMPRM can achieve SLAM with the best performance. Therefore, 3DMNDT
has the potential for the SLAM application.

##  VI Conclusions

Under the framework of NDT, this paper proposes a novel method for the multi-
view registration. To the best of our knowledge, this is the ﬁrst time the NDT
is applied to multi-view registration. Instead of space division, it utilizes
k-means clustering algorithm to cluster data points and then computes NDT for
each cluster. Compared with regular division, the K-means clustering can avoid
unbalanced division of data points. Then, it formulates the multi-view
registration problem by the NDT-based likelihood function, which is maximized
by the Lie algebra solver. Different from the Newton’s solver, the Lie algebra
solver does not require to calculate Jacobian matrix and Hessian matrix, which
is tedious and error-prone in 3D case. Experimental results tested on bench
mark data sets illustrate that the proposed method can achieve state-of-the-
art performance. What’s more, it can be applied to SLAM.

## Acknowledgment

This work is supported by the National Natural Science Foundation of China
under Grant No. 61573273.

## References

  * [1] R. Sandhu, S. Dambreville, and A. Tannenbaum, “Point set registration via particle filtering and stochastic dynamics,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 32, no. 8, pp. 1459–1473, 2009\.
  * [2] J. Yang, H. Li, D. Campbell, and Y. Jia, “Go-icp: A globally optimal solution to 3d icp point-set registration,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 38, no. 11, pp. 2241–2254, 2015.
  * [3] J. Ma, J. Zhao, J. Jiang, H. Zhou, and X. Guo, “Locality preserving matching,” _International Journal of Computer Vision_ , vol. 127, no. 5, pp. 512–531, 2019.
  * [4] Z. Jiang, J. Zhu, Y. Li, J. Wang, Z. Li, and H. Lu, “Simultaneous merging multiple grid maps using the robust motion averaging,” _Journal of Intelligent & Robotic Systems_, vol. 94, no. 3-4, pp. 655–668, 2019.
  * [5] W. Tabib, C. O’Meadhra, and N. Michael, “On-manifold gmm registration,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 4, pp. 3805–3812, 2018\.
  * [6] Z. Jiang, J. Zhu, Z. Lin, Z. Li, and R. Guo, “3d mapping of outdoor environments by scan matching and motion averaging,” _Neurocomputing_ , vol. 372, pp. 17–32, 2020.
  * [7] A. Rasoulian, R. Rohling, and P. Abolmaesumi, “Group-wise registration of point sets for statistical shape models,” _IEEE transactions on medical imaging_ , vol. 31, no. 11, pp. 2025–2034, 2012.
  * [8] N. Ravikumar, A. Gooya, S. Çimen, A. F. Frangi, and Z. A. Taylor, “Group-wise similarity registration of point sets using student’s t-mixture model for statistical shape models,” _Medical image analysis_ , vol. 44, pp. 156–176, 2018.
  * [9] P. J. Besl and N. D. McKay, “Method for registration of 3-d shapes,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 14, no. 2, pp. 239–256, 1992.
  * [10] P. Biber and W. Straßer, “The normal distributions transform: A new approach to laser scan matching,” in _Proceedings 2003 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2003)(Cat. No. 03CH37453)_ , vol. 3. IEEE, 2003, pp. 2743–2748.
  * [11] M. Magnusson, “The three-dimensional normal-distributions transform: an efficient representation for registration, surface analysis, and loop detection,” Ph.D. dissertation, Örebro universitet, 2009.
  * [12] S. Rusinkiewicz and M. Levoy, “Efficient variants of the icp algorithm.” in _3dim_ , vol. 1, 2001, pp. 145–152.
  * [13] D. Chetverikov, D. Stepanov, and P. Krsek, “Robust euclidean alignment of 3d point sets: the trimmed iterative closest point algorithm,” _Image and vision computing_ , vol. 23, no. 3, pp. 299–309, 2005.
  * [14] A. Myronenko and X. Song, “Point set registration: Coherent point drift,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 32, no. 12, pp. 2262–2275, 2010.
  * [15] B. Jian and B. C. Vemuri, “Robust point set registration using gaussian mixture models,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 33, no. 8, pp. 1633–1645, 2010.
  * [16] W. Gao and R. Tedrake, “Filterreg: Robust and efficient probabilistic point-set registration using gaussian filter and twist parameterization,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 11 095–11 104.
  * [17] I. Sipiran and B. Bustos, “Harris 3d: a robust extension of the harris operator for interest point detection on 3d meshes,” _The Visual Computer_ , vol. 27, no. 11, p. 963, 2011.
  * [18] R. B. Rusu, N. Blodow, and M. Beetz, “Fast point feature histograms (fpfh) for 3d registration,” in _2009 IEEE international conference on robotics and automation_. IEEE, 2009, pp. 3212–3217.
  * [19] Y. Guo, F. Sohel, M. Bennamoun, M. Lu, and J. Wan, “Rotational projection statistics for 3d local surface description and object recognition,” _International journal of computer vision_ , vol. 105, no. 1, pp. 63–86, 2013\.
  * [20] H. Lei, G. Jiang, and L. Quan, “Fast descriptors and correspondence propagation for robust global point cloud registration,” _IEEE Transactions on Image Processing_ , vol. 26, no. 8, pp. 3614–3623, 2017.
  * [21] Y. Aoki, H. Goforth, R. A. Srivatsan, and S. Lucey, “Pointnetlk: Robust & efficient point cloud registration using pointnet,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 7163–7172.
  * [22] Y. Wang and J. M. Solomon, “Deep closest point: Learning representations for point cloud registration,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2019, pp. 3523–3532.
  * [23] C. Choy, W. Dong, and V. Koltun, “Deep global registration,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 2514–2523.
  * [24] Y. Chen and G. Medioni, “Object modelling by registration of multiple range images,” _Image and vision computing_ , vol. 10, no. 3, pp. 145–155, 1992\.
  * [25] R. Bergevin, M. Soucy, H. Gagnon, and D. Laurendeau, “Towards a general multi-view registration technique,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 18, no. 5, pp. 540–547, 1996.
  * [26] J. Zhu, Z. Li, S. Du, L. Ma, and T. Zhang, “Surface reconstruction via efficient and accurate registration of multiview range scans,” _Optical Engineering_ , vol. 53, no. 10, p. 102104, 2014.
  * [27] Y. Tang and J. Feng, “Hierarchical multiview rigid registration,” in _Computer Graphics Forum_ , vol. 34, no. 5. Wiley Online Library, 2015, pp. 77–87.
  * [28] S. Krishnan, P. Y. Lee, J. B. Moore, S. Venkatasubramanian _et al._ , “Global registration of multiple 3d point sets via optimization-on-a-manifold.” in _Symposium on Geometry Processing_ , 2005, pp. 187–196.
  * [29] X. Mateo, X. Orriols, and X. Binefa, “Bayesian perspective for the registration of multiple 3d views,” _Computer Vision and Image Understanding_ , vol. 118, pp. 84–96, 2014.
  * [30] V. M. Govindu and A. Pooja, “On averaging multiview relations for 3d scan registration,” _IEEE Transactions on Image Processing_ , vol. 23, no. 3, pp. 1289–1302, 2013.
  * [31] J. Sola, J. Deray, and D. Atchuthan, “A micro lie theory for state estimation in robotics,” _arXiv preprint arXiv:1812.01537_ , 2018.
  * [32] F. Arrigoni, B. Rossi, and A. Fusiello, “Global registration of 3d point sets via lrs decomposition,” in _European Conference on Computer Vision_. Springer, 2016, pp. 489–504.
  * [33] F. Arrigoni, B. Rossi, P. Fragneto, and A. Fusiello, “Robust synchronization in so (3) and se (3) via low-rank and sparse matrix decomposition,” _Computer Vision and Image Understanding_ , vol. 174, pp. 95–113, 2018.
  * [34] R. Guo, J. Zhu, Y. Li, D. Chen, Z. Li, and Y. Zhang, “Weighted motion averaging for the registration of multi-view range scans,” _Multimedia Tools and Applications_ , vol. 77, no. 9, pp. 10 651–10 668, 2018.
  * [35] C. Jin, J. Zhu, Y. Li, S. Pang, L. Chen, and J. Wang, “Multi-view registration based on weighted lrs matrix decomposition of motions,” _IET Computer Vision_ , vol. 13, no. 4, pp. 376–384, 2018.
  * [36] G. D. Evangelidis and R. Horaud, “Joint alignment of multiple point sets with batch and incremental expectation-maximization,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 40, no. 6, pp. 1397–1410, 2017\.
  * [37] J. Zhu, Z. Jiang, G. D. Evangelidis, C. Zhang, S. Pang, and Z. Li, “Efficient registration of multi-view point sets by k-means clustering,” _Information Sciences_ , vol. 488, pp. 205–218, 2019.
  * [38] J. Zhu, R. Guo, Z. Li, J. Zhang, and S. Pang, “Registration of multi-view point sets under the perspective of expectation-maximization,” _IEEE Transactions on Image Processing_ , vol. 29, pp. 9176–9189, 2020.
  * [39] A. Das and S. L. Waslander, “Scan registration with multi-scale k-means normal distributions transform,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2012, pp. 2705–2710.
  * [40] A. Nuchter, K. Lingemann, and J. Hertzberg, “Cached kd tree search for icp algorithms,” in _Sixth International Conference on 3-D Digital Imaging and Modeling (3DIM 2007)_. IEEE, 2007, pp. 419–426.
  * [41] C. Forster, L. Carlone, F. Dellaert, and D. Scaramuzza, “On-manifold preintegration for real-time visual–inertial odometry,” _IEEE Transactions on Robotics_ , vol. 33, no. 1, pp. 1–21, 2016.
  * [42] J. Zhu, L. Zhu, Z. Li, C. Li, and J. Cui, “Automatic multi-view registration of unordered range scans without feature extraction,” _Neurocomputing_ , vol. 171, pp. 1444–1453, 2016.