Toward a Unified Framework for Point Set Registration
Feiran Li, Kent Fujiwara and Yasuyuki Matsushita

Abstract— Point set registration plays a critical role inrobotics and computer vision. Early methods considered registration as a purely geometric problem, presenting excellent
extensibility for various tasks due to their explicit handling ofcorrespondences; statistical methods were later introduced tohandle noise. However, the two categories of algorithms have
evolved independently without sharing much in common. In
this paper, we leverage the concept of information geometry totheoretically unify the two classes together by interpreting themas the same operation but in different spaces associated withrespective metrics. Moreover, based on the proposed unification,
we also develop a novel bandwidth estimation strategy to solve
the long-standing problem of statistical registration algorithms,
and demonstrate its theoretical and practical advantages overdeterministic annealing, the most commonly used technique. We
also present a case study to show how geometric and statistical
approaches can benefit from each other.

I. I NTRODUCTION
Pairwise point set registration is a central problem incomputer vision and graphics, which plays an importantrole in 3D reconstruction, shape retrieval, medical imaging,
and robot vision. It is formulated as finding the optimal
rigid transformation T, consisting of rotation R and translation t, between two point sets A = {a1 , a2 , . . . , am } andB = {b1 , b2 , . . . , bn }.
The difficulty of registration lies in the unknown correspondences between bj and ai . As pioneering approaches,
geometric algorithms such as Iterative Closest Point (ICP) [1]first proposed to tackle this problem by alternating between establishing correspondences and optimizing T. Other
methods based on statistics [2], [3], graph matching [4],feature correspondence [5], and deep learning [6] were later
introduced to address various challenges. In this work, we
focus on the geometric and statistical approaches since they
are used by a vast majority of registration tasks.
While statistical methods excel at coping with various typesof noise, geometric ones provide more choices in various
conditions due to their explicit handling of correspondences.
Despite the maturity of these methods, there still exist somelong-standing problems. First, the two kinds of algorithms
evolved independently. As a result, the relationship between
the two classes is unclear, preventing them from sharing
techniques for mutual benefit. Second, existing statisticalregistration methods lack effective strategies to estimate
the bandwidths, and heavily depend on empirical tuning∗ This research was supported by NII CRIS collaborative research program
operated by NII CRIS and LINE Corporation.
∗ Feiran Li and Yasuyuki Matsushita are with Osaka University. Kent
Fujiwara is with LINE Corporation.
∗ {li.feiran, yasumat}@ist.osaka-u.ac.jp, kent.fujiwara@linecorp.com

without any theoretical justifications [7]. Consequently, theirperformance is hindered from further enhancement.
In this work, we aim at developing a unified framework as
well as a reasonable bandwidth estimation strategy to tackle
the aforementioned problems. For unification, we show thatthe relation between geometric and statistical registration
algorithms can be explained from an information geometricperspective, where both geometric and statistical approaches
are regarded as conducting the same operation but in different
spaces. For bandwidth estimation, we show that the objectivefunction developed for registration intrinsically serves as a
guidance for estimating bandwidths; hence we do not needto tune them empirically. We summarize our contributions asfollows:• A theoretical unification that bridges the geometric andstatistical registration approaches for mutual benefit.
• A bandwidth estimation strategy with theoretical andpractical advantages over existing approaches.
It is noteworthy that there already exist some attempts
to bridge the geometric approaches with the statistical ones
via expectation-maximization (EM) [7], [8]. However, thefact that the underlying metric of the EM algorithm is the
Kullback-Leibler divergence (KL-divergence) prevents theseworks from explaining other divergence-minimization-based
methods [9], [10]. As a result, such a unification cannot
generalize to all the cases. On the other hand, our explanation
suits all cases well without exception.
II. R ELATED WORK
Geometric and statistical methods possess their respective
merits in solving the point set registration problem. Specifically, owing to the explicit handling of correspondences,geometric methods are more flexible in tasks such as
symmetric registration and outlier rejection. On the other
hand, statistical ones are more accurate under conditions suchas point-wise varying noises or clustering-based registration.
We summarize their common objectives in Table I, where
MM and KDE respectively stand for mixture model (a.k.a.,
mixture family distribution) and kernel density estimation.
a) Geometric approaches for registration: By formulating the objective as the orthogonal Procrustes problem [11]with unknown correspondences, ICP [1] leads the development of geometric approaches for point set registration.
Some variants of it present alternatives on establishingcorrespondences. For example, EM-ICP [7] relaxes the oneto-one correspondence rule to soft assignment to deal withinaccurate initialization. Gold et al. [12] convert the soft
assignment matrix into a doubly stochastic one to allow
symmetric registration. Trimmed-ICP [13] robustifies theTABLE I: Summary of geometric and statistical registration
algorithms.

Remark 1: The procedure of maximizing likelihood isequivalent to minimizing KL-divergence: argmax P (x; θ) =θGeometric
Statistical
Statistical

Problem formulation
Orthogonal Procrustes problem
w/ unknown correspondenceDistribution-fittingof MM
Divergence-minimization
between MM/KDE

Exampleargmin KL (P (x; θ? ) |P (x; θ)); where θ is the estimated

ICP [1]parameter and θ? is the ground-truth one.
Per Remark 1, the objective function of distribution-fitting
methods can be rewritten as minimizing divergence:YargmaxP (bj ; A, T, θ) ≡T
(1)argmin KL (P (B; θb ) |P (A; T, θa )) ,CPD [2]
GMMReg [3]standard ICP algorithm by trimming out probable miscorrespondences. There also exist other extensions that study
global optimization and efficiency. For instance, Yang etal. [14] propose a Branch-and-Bound search strategy toguarantee global optimization. Greenspan and Yurick [15]use the approximated k-d tree for reducing the computational
cost of matching. The point-to-plane [16], [17] and planeto-plane [18] variants substitute the point-to-point correspondence for faster convergence and symmetric registration.
b) Statistical approaches for registration: Statisticalmethods treat point sets as mixture family distributions tocope with noise. Some works consider it as a distributionfitting problem. For example, Myronenko et al. [2] consider
one point set as a Gaussian mixture model (GMM) andanother one as samples drawn from it. Consequently, it casts
the registration problem to estimating the parameters of the
distribution, which is solvable via the EM algorithm. Anotherthread of work fits MMs to both point sets, and the optimaltransformation can be estimated by minimizing the divergence
between them. For example, Jian and Vemuri [8] recommend
minimizing the L2-divergence between GMMs. Other distributions and divergences are also explored, such as Student’s tdistribution [19], hybrid mixture model [20], Jensen-Shannon
divergence [10] and Cauchy-Schwarz divergence [9].
c) Attempts to unify geometric and statistical approaches: There already exists literature that attempts tounify geometric and statistical approaches. For example,
Granger and Pennec [7] interpret ICP as minimizing the sumof Mahalanobis distances between Gaussian distributions.
Jian and Vemuri [8] explain ICP as an approximation ofminimizing KL-divergence between two GMMs. Segal etal. [18] use maximum likelihood estimation (MLE) tomaximize the sum of probabilities of distances between paired
points, which actually leads to the same conclusion drawn
by [8] but from a different perspective. The main difference
between these existing works and our proposal is that,while they treat ICP as a particular case of their respective
objective functions without referring to its variants, we
provide thorough theoretical explanations to unify geometricand statistical methods without any exception.
III. A UNIFIED FRAMEWORK FOR GEOMETRIC ANDSTATISTICAL APPROACHES

We first need to unify distribution-fitting and divergenceminimization methods within the statistical class before
unifying the geometric and the statistical ones.

θ

T

P
where
P (A; T, θa ) =αi P (ai ; T, θai ) and P (B; θb ) =P
βj P (bj ; θbj ) are MMs constructed from P
point sets A and
B. We hereafter
denote
them
as
P
(A)
=
αi P (ai ) andP
P (B) =
βj P (bj ) for notational simplicity. Now, all thestatistical methods can be seen as a divergence-minimization
problem, which lays the foundation for our unification.
KL-divergence is not the only choice for formulating theobjective function. In fact, as employed by numerous related
works [8]–[10], there are many off-the-shelf divergences
for selection. However, not all of them are suitable for theregistration problem. As a naive rule, we desire the selected
divergence to be invariant when the points are encoded indifferent ways. Specifically, given a mapping h, we want thedivergence to satisfy
D (P |Q) = D (h (P ) |h (Q)) .

(2)

Eq. (2) directly leads to the theory of f -divergence, which contains a group of famous divergences such as KL-divergence,
total variation, and all of the α-divergences. We present its
formal definition in Lemma 1.
 RLemma 1: The f -divergence: Df (P |Q) , Ω pf pq dµ
is an invariant measurement, where µ is the referencedQ
distribution satisfying p = dPdµ and q = dµ , and f is asmooth convex function. It is convex w.r.t. both P and Q.
A detailed introduction together with proof of Lemma 1
can be found in P54-P65 of [21]. Although the f -divergenceis a reasonable choice to formulate the objective function,
it is seldom analytically computable when the incorporated
distributions are MMs. Therefore, we need to find a solvable
upper bound for optimization purposes.
PM
Given two MMs written in the form of P = i αi Pi
PN
and Q = j βj Qj , where Pi and Qj are the respective
component distributions of P and Q, and αi and βj are theirweights; we can treat them as having the same number ofcomponentsPM N
 P = i αi,n Pi,n,PM N
Q = j βj,m Qj,mαi,1 = · · · = αi,n = αNi





(3)


βj

 βj,1 = · · · = βj,m = M
where
.



Pi,1 = · · · = Pi,n = Pi






Qj,1 = · · · = Qj,m = Qj

Euclidean distance

Algorithm 1: The unified algorithm for point set
registration
Input: Point sets A and B, initial transformation T
While not converged:
3
3• Step 1: Find a mapping ε : R → M (or R ) for
each point ai and bj under current T.
3• Step 2: Establish correspondences on M (or R )under current T.
• Step 3: Minimize the sum of divergences (orEuclidean distances) w.r.t. T.

(a) Geometric

Owing to its convexity, we can define an upper bound of the
f -divergence Df (P |Q) in the summation form via Jensen’s
inequality [22]


M,N
M,N
X
X
Df (P |Q) = Df 
αi,n Pi,n |βj,m Qj,m i,n6

M,N,M,N
X

j,m

(4)

Df (αi,n Pi,n |βj,m Qj,m ) .

i,n,j,mMoreover, since MMs are invariant to permutation of indices,
we can always tighten the upper bounds in the form of
Df (P |Q) 6

M,N,M,NX

Df αi,n Pi,n |βs(i) Qs(i)

i,n,j,m

6

M,N,M,NX

πij Df (αi,n Pi,n |βj,m Qj,m )(5)

i,n,j,m=

M,N,M,NX
i,n,j,m
πij Dfβjαi
Pi,n | Qj,m ,NM

where βs(i) Qs(i) denotes the nearest component
P of Q w.r.t.
αi,n Pj,m , and πij is a weight coefficient s.t. j πij = 1.
Since the MMs generated by the two point sets commonly
1obey αi = βj = M, by applying Eq. (5) to the registration
problem, the objective function can be written as
M,N,M,NX
1
1J=
πij DfP (ai ) |P (bj ) . (6)MN
NMi,n,j,mMoreover, according to the invariance property of the f divergence, we can equally scale each component distribution
by M N without changing the objective function. Consequently, the above equation is simplified to
J=

M,N
X

πij Df (P (ai ) |P (bj )) .

(7)

i,j

Recalling that the objective functions of ICP in its generalized form (i.e., soft assignment) is written asX2J=
πij kT ◦ ai − bj k2 ,(8)
we can see that it is very similar to Eq. (7), except that
the Euclidean distance is replaced by the f -divergence. To

Fisher Information Metric(b) StatisticalFig. 1: Conceptual illustration of how geometric and statisticalmethods work. While geometric approaches operate on
Euclidean space, statistical ones work on statistical manifold.

have a more comprehensive understanding of what statisticalmethods actually do, we introduce another property of the
f -divergence shown in Lemma 2.
Lemma 2: The f -divergence naturally introduces geometric structures to P and Q, where fisher information metric(FIM) is the associated Riemannian metric. Also, we canalways use a convex function f satisfying f (1) = 0, f 0 (1) =0 and f 00 (1) = 1 to formulate a standard f -divergence.
Lemma 2 in fact leads to the theory of informationgeometry [21]. It studies statistics with geometric tools bytreating distributions as points on a Riemannian manifold M.
For proof:
Proof: Given two neighboring points P (x; θ) andP (x; θ + δθ) (hereafter denoted as P (θ) and P (θ + δθ) for
notational simplicity) on M, we can obtain the local formof the f -divergence:Z
P (θ + δθ)
Df (P (θ) |P (θ + δθ)) = P (θ) f
. (9)
P (θ)
By expanding the term f (·) with Taylor expansion to the
2nd order and using the properties f (1) = 0, f 0 (1) = 0 and
f 00 (1) = 1, Eq. (9) can be rewritten into
 


Z
P (θ)
∂f
1 T ∂2f
Df = P (θ) f+δθ + δθ
δθ
P (θ)
∂θ
2
∂θ2


Z
∂2f
1
P (θ) · 0 + 0 + δθT 2 δθ
=
2
∂θ
Z 
 


1
1
∂
∂
= δθT
·
·
P (θ)
P (θ)
δθ
2
P (θ)
∂θ
∂θ
Z

1
∂ log P (θ) ∂ log P (θ)
= δθT
P (θ)
δθ
2
∂θ
∂θ
1
= gij dθi dθj ,
2
(10)
where
Einstein
Summation
form
is
used
and
g
=
ij
R
P (θ) ∂ log P (θ)
P (θ) ∂ log
isexactly
the
FIM.
∂θi
∂θj
Since the last line of Eq. (10) is a general notation for
distance, we can also express the euclidean distance in thesame form by changing gij to the Kronecker-delta. Therefore,
we can now express the objective functions of both statisticaland geometric methods in the unified form ofπk k i jJ=g dθ dθ ,(11)
2 ij k kwhere k denotes the indices of pairs of points; and gij encodes
distinct metrics for statistical and geometric methods.
With this unified objective function, we can conclude thatthe two classes of algorithms, in fact, conduct the same
operations in different spaces associated with different metrics.
Specifically, as illustrated in Fig. 1, they all actually do ICP(or its variants) but in distinct spaces.
One merit of the unified explanation is that, we are freed
from having to select between geometric and statistical registration algorithms. I.e., we can always develop an algorithm
that imitates ICP and its variants on handling correspondences
while keeping the algorithm statistical to handle noise. We
summarize the unified algorithm in Algorithm 1, where ε isan endomorphism (R3 → R3 ) for geometric methods and ahomomorphism (R3 → M) for statistical counterparts.
IV. A N INTRINSICALLY CONTAINED METHOD FOR
BANDWIDTH ESTIMATION STRATEGY

We now put Algorithm 1 into practice. A major problem
faced by many divergence-minimization-based approaches
lies in how to define an effective mapping ε : R3 → M.
Fortunately, the aforementioned f -divergence intrinsicallycontains an approach for doing so. To maintain consistency
with other registration literature, and without loss of generality,
we hereafter assume that all the points within a given point
set are contaminated by the same isotropic Gaussian noise.
Since the expectation of each Gaussian distribution can benaturally determined as the point position in the Euclidean
space, the only remaining problem related to the mapping
pertains to the covariances Σa and Σb . In the case ofclustering-based approaches [23], [24], they are easy toestimate and can be treated as known parameters. However,in general cases, the mapping is considered as a bijection,
which leads to the bandwidth estimation problem of KDE.
Empirically, the bandwidth of a general KDE can be
estimated beforehand and fixed with rule-of-thumb methods [25]. However, in the context of the registration problem,a dynamically updated one would be more suitable since therelative position of the two point sets varies according to thetransformation T.
The heuristic method of deterministic annealing is a popular
tool for setting the bandwidths Σa and Σb [3], [8], [9].
Its philosophy is based on the asymptotic behaviors as
studied by Chui et al. [26]. In detail, Σa and Σb control thecorrespondences πij between point sets. The correspondences
are vague when Σa and Σb are large, and become asymptotically one-to-one when they approach 0. Therefore, in orderto take advantage of both soft assignment and one-to-one
correspondences to avoid local minima and obtain accurate
registration, deterministic annealing gradually decreases Σaand Σb from large initial values until convergence. Althoughpromising performances have been observed, we hesitate to
do so since bandwidth estimation is an independent problem
in statistics with more effective solutions.

We claim that the objective function developed for registration intrinsically acts as a guidance for dynamically setting
the bandwidths. Specifically, it reminds us of the bandwidth
estimation technique based on risk functions [27].
Remark 2: A risk function L is defined as the expectationof a loss f : L = E [f ]. Applying to bandwidth estimation of
KDE, the optimal bandwidth of a KDE p can be retrieved byminimizing E [f (p, q)], where q is the underlying groundtruth distribution for reference.
According to the definition in Remark 2, the f -divergencein fact naturally
one risk function in the formh formulates
iof Df = Ep f pq . Therefore, if we alternately fix thepoint P (ai ) or P (bj ) as the reference, the f -divergence also
serves as an objective function for estimating the bandwidthof the other one.
If the reference distribution q is the ground-truth, a plausible
bandwidth of the variable distribution p would be the onethat minimizes the difference between them. However, for theregistration problem, since q is just temporarily fixed without
any guarantees to be faithful, directly tuning p towards it
would unavoidably introduce large variance for smoothing
purposes. This problem is known as the bias-variance-tradeoff,
and we employ the minimizing entropy criterion proposedby Jiang et al. [28] as a solution. Specifically, we can limit
the variance of a KDE by adding in a regularization term of
its entropy, which leads to our final objective function:X
J (T, Σa , Σb ) =πij Df (P (ai ) |P (bj ))(12)+ λ1 H (P (A)) + λ2 H (P (B)) ,where H (·) is the entropy, and λ1 and λ2 are the regularization parameters. As commonly used in other registration
algorithms, Eq. (12) can be facilely minimized via alternating
optimization, which update T, Σ and π alternately.
V. H OW UNIFICATION HELPS MUTUAL BENEFIT:
A CASE STUDY
Based on the aforementioned unification, we can always
design an algorithm that takes advantage of both statistical andgeometric methods, namely, an algorithm based on divergenceminimization with explicitly maintained correspondences. In
this section, we take the symmetric registration task as an
example to demonstrate its effectiveness.
A. Problem formulation
Assuming that the two point sets consist of a similar amountof points and no preliminary information of the registration
order (A → B or B → A) is given, symmetric registration
aims at alleviating the effects of registration order on the
final results. That is, we want A → B and B → A to returnnearly the same estimations.
This problem is easy to solve for geometric methods. For
example, as used in [12], [29], [30], it is a common practice
to convert the aforementioned soft assignment matrix π toan (asymptotic) doubly stochastic one in each iteration. In
detail, the property of doubly stochastic matrices that eachrow and column sums to 1 helps to ignore the normalization

100

ICP
RPM
CPD
SVR
Ours

25
20
15
10

T-Err (1e-2)

T-Err (1e-2)

30

5
0

0.1

0.2

0.3

0.4

0.5

Covariance of noise (1e-2)

ICP
RPM
CPD
SVR
Ours

80
60
40
20
0

0.1

0.2

0.3

0.4

Ratio of outliers

0.5

Fig. 2: An example of the synthetic data contaminated by
Gaussian noise and outliers. Left: point clouds in their initialposes. Right: Our registration result.

Fig. 3: Comparison w.r.t. varying levels of noise and ratios
of outliers. Each experiment consist of 10 different setups.

TABLE II: Errors on synthetic data (both rotational andtransformational errors are in 1e − 2 scale).

We use the symmetrized KL-divergence to specify the f divergence for maintaining symmetry. Also, in order to showthat the tools developed for geometric methods can be directly
applied to statistical ones, we follow RPM [7] and use softmax
for soft assignment and Sinkhorn iteration [35] for converting
the assignment matrix π to doubly stochastic. For each pair ofpoint sets in both synthetic and realistic datasets, we conduct
the registration from both A → B and B → A. We use
R−Err = I − RTgt Resti F , t−Err = ktesti − tgt k2 andT−Err = I − T−1esti Tgt F as the respective error metrics
for rotation, translation and transformation. For parameter
setup, the regularization coefficients for bias-variance-tradeoff
in our method are set to λ1 = λ2 = 50. The annealing
parameter is set to 0.95 for RPM and 0.9 for SVR, withwhich we observe better performances. For CPD, the weightof uniform distribution is set to 0.1. The maximum iteration is
set to 100 for all the algorithms except SVR, which is set to
10 for its slow clustering procedure. In fact, we observe thatonce the SVR algorithm succeeds, the number of iterations
shows minor effects on its final accuracy.
a) Tests on synthetic data: We first test how different
algorithms perform under noisy conditions. For experimental
setup, we base the tests on the bunny and the dragon pointsets given in [36], [37], and a teddy from Free3D2 . They
are all down-sampled to around 1000 points for a reasonableruntime, as shown in Fig. 2. We add isotropic Gaussian noises
with standard deviation σ = 0.0015I and σ = 0.0025I, as
well as 10% of outliers to the two point sets respectively.
The registration is repeated 20 times on each pair of pointsets. For each trial, a random ground-truth rotation R drawn
from Euler angles in the range of (0◦ , 60◦ ) along each axis,
and a random translation is used. The results are shown in
Table II. As we can see, our statistical symmetric registration
method can champion both the rotational and transformational
estimations, although it performs negligibly inferior to CPD
on estimating the translations.
We also test the robustness of our algorithm w.r.t. different
ratios of outliers and levels of noises. For setup, we use
the aforementioned bunny point set and fix the ground-truth
rotation to 30◦ in Euler angles. For outliers, we use the 10
ratios uniformly distributed within 5% to 50%. For noise,
we vary the covariances of Gaussian noise from 1e − 3I to
5e − 3I with a step size of 5e − 4I. We repeat the experiment

Mean
ICP
RPM
CPD
SVR
Ours

R err
39.49
4.92
1.67
11.73
0.93

t err
6.93e-4
7.04e-43.30e-4
6.02e-4
3.73e-4Median
T err
39.49
4.93
1.67
11.73
0.93

R err
2.91
4.29
1.32
4.21
0.96

t err
6.71e-4
6.96e-4
3.24e-44.08e-4
3.56e-4

T err
2.91
4.29
1.32
4.21
0.96

axis of π, hence the registration order is perfectly balanced.
On the other hand, dealing with this task is a weakness ofexisting statistical methods. In detail, the “A as model and
B as data” principle of distribution-fitting-based methods [2],[20] naturally requires a pre-defined registration order. For
the divergence-minimization-based methods, although we canrely on some symmetric divergences [9], [10], there still lacks
a unified solution when using general divergences.
Our unification helps the statistical algorithms to strengthen
symmetric registration. In detail, if we treat each component
distribution as a point on M, we can naturally follow thegeometric methods and convert the assignment matrix π todoubly stochastic in each iteration. The only difference isthat, while geometric methods use the Euclidean distance,
the statistical ones must replace it with a divergence.
B. Experiments
We select 4 other candidates for comparison: ICP [1],
RPM [12], CPD [2] and SVR [31]. Among them, ICP solves
the registration problem in a pure geometric manner. RPM
converts the correspondences of ICP to doubly stochastic torealize symmetric registration. CPD and SVR respectively
represent the distribution-fitting and divergence-minimizationbased statistical approaches. Regarding the implementations
used in the experiments, ICP and RPM are implemented
by ourselves, CPD is provided by the authors, and SVR isfrom the ProbReg Python library1 . For all of the following
experiments, we assume the datasets are not under extreme
conditions (i.e., severe outliers, low overlaps, or significant distinction between numbers of points) since they are commonly
considered as independent tasks addressed by extending the
4 base approaches with other techniques [13], [32]–[34].
1 ProbReg: http://probreg.readthedocs.io/. Last accessed on
March 24, 2021.

2 Free3D http://www.free3d.com. Retrieved March 24, 2021.

12

2.0

DA
Ours

1.5
1.0
0.5
0.95

TABLE III: Errors on the laser scan data.
Mean
ICP
RPM
CPD
SVR
Ours

R err
2.11e-2
5.06e-2
2.08e-2
6.69e-2
1.66e-2

t err
0.61
0.11
0.15
0.48
0.12

10
8
6
4
2

0.96

0.97

0.98

0.99

Annealing coefficient s

Fig. 4: An example of the ETH Hauptgebaude dataset. Somepoints are cropped for better visualization. Left: the pointsets in their initial poses. Right: our registration result.

Trials countT-Err (1e-2)

2.5

0

0.2 0.4 0.6 0.8 1.0 1.2 1.4

Transformational err (1e-2)

Fig. 5: Comparison of deterministic annealing and ourbandwidth estimation strategy. Left: Accuracy w.r.t. different
annealing parameters. The orange line is achieved by fixing
the tradeoff parameters to λ1 = λ2 = 50. Right: Error
distribution with λ1 = λ2 uniformly drawn from [1, 60].

Median
T err
0.61
0.11
0.15
0.49
0.12

R err
1.90e-2
4.65e-2
1.93e-2
4.21e-2
1.56e-2

t err
0.55
0.11
0.15
0.36
0.12

T err
0.61
0.11
0.15
0.49
0.12

10 times with each ratio and each level. The mean errors areplotted in Fig. 3. As shown, for the case of different levels ofnoise, our proposal in general performs better than CPD. For
the case of distinct ratios of outliers, our method can stablypresent higher accuracy compared with the others.
b) Tests on real data: We also conduct experiments onlaser scan dataset to study the real-world performances. In
detail, we base the tests on the first 10 scans of the ETH
Hauptgebaude dataset [38], as shown in Fig. 4. This dataset ischallenging as it consists of repetitive elements. Since some ofthe original ground-truth rotations are either extremely largethat none of the candidates can give reasonable estimations,
we manually reset the rotations to the range of (0◦ , 60◦ ) in
Euler angles along each axis. We also change the weightof the uniform distribution of CPD to 0 as we found it
barely succeeds with other numbers. Again, each point set
is down-sampled to around 1000. The results are shownin Fig. III. As can be seen, our method outperforms all
the others on estimating rotations. For translational andtransformational errors, although RPM presents slightly better
results, the relative differences between them and ours are
small compared to the improvements on rotations.
VI. P ERFORMANCE OF THE PROPOSED BANDWIDTHESTIMATION STRATEGY

In this section, we study the performance of our proposed
bandwidth estimation strategy. For experimental setup, we use
the aforementioned synthetic bunny with a fixed ground-truth
rotation of (30◦ , 30◦ , 30◦ ) in Euler angles; and keep noise,
outliers, and parameter settings the same as mentioned above.
a) Comparisons between deterministic annealing andour proposal: We study the relative performance between
the commonly used deterministic annealing technique andour proposal to show the superiority of our bandwidth
estimation strategy. Specifically, given initial bandwidths,

the deterministic annealing technique multiplies it with ashrinkage coefficient s after each iteration. In our implementation, we observe that the registration performances become
significantly unstable when s is within the range of (0.9, 0.95),
and it hardly succeeds when s ∈ (0, 0.9). Therefore, we
limit s to be drawn from (0.95, 1) to make the comparison
meaningful. The results are reported in Fig. 5. As shown in
the figure, the performance of deterministic annealing is notstable as minor changes on s led to conspicuous shifts in theestimated results. Moreover, its accuracy is still significantly
inferior to our strategy.
b) Stability w.r.t. different λ: As the only free parametersin our proposal, it is meaningful to study how different λ1and λ2 for bias-variance-tradeoff can affect the performance.
We again simply assume λ1 and λ2 to share the same value:
λ1 = λ2 = λ. For testing purposes, we uniformly pick 20
trials with λ ∈ [1, 60]. The final errors are summarized in thehistogram in Fig. 5. As can be seen, our bandwidth estimation
strategy is stable w.r.t. varying λ within a broad tuning range.
In detail, the apparent dominant error suggests that a large
scope of different λ can eventually lead the algorithm to thesame estimation. Moreover, compared to the performance ofthe deterministic annealing mentioned above, the difference
between the best and worst estimations of our proposal is
only in 10−3 magnitude.
VII. D ISCUSSION AND CONCLUSION
In this work, we present a unified framework by explaining
geometric and statistical point set registration algorithms
as conducting the same operations in different spaces. We
also derive a bandwidth estimation strategy for generaldivergence-minimization-based methods, which demonstrates
its effectiveness over the existing heuristic method. Also, by
taking the symmetric registration task as an example, weshow how to combine the strengths of both classes, whose
effectiveness is demonstrated in various experiments.
For future work, we plan to combine the statisticalregistration methods with other geometric ones to deal withcases such as partial overlaps and sparse-to-dense registration.
Another direction lies in exploring how to select a suitablef -divergence for the registration problem, which remains an
open problem in information science.

R EFERENCES[1] P. J. Besl and N. D. McKay, “Method for registration of 3-d shapes,”in Transactions on Pattern Analysis and Machine Intelligence (PAMI),
vol. 14, pp. 239–256, 1992.
[2] A. Myronenko and X. Song, “Point set registration: Coherent point drift,”
Transactions on Pattern Analysis and Machine Intelligence (PAMI),
vol. 32, no. 12, pp. 2262–2275, 2010.
[3] B. Jian and B. C. Vemuri, “Robust point set registration using gaussian
mixture models,” Transactions on Pattern Analysis and MachineIntelligence (PAMI), vol. 33, no. 8, pp. 1633–1645, 2010.
[4] F. Zhou and F. De la Torre, “Factorized graph matching,” in Proceedingsof Conference on Computer Vision and Pattern Recognition (CVPR),
pp. 127–134, 2012.
[5] S. Gumhold, X. Wang, and R. S. MacLeod, “Feature extraction frompoint clouds,” in Proceedings of International Meshing Roundtable,
Sandia National Laboratories, 2001.
[6] G. Elbaz, T. Avraham, and A. Fischer, “3d point cloud registration for
localization using a deep neural network auto-encoder,” in Proceedingsof Conference on Computer Vision and Pattern Recognition (CVPR),
pp. 4631–4640, 2017.
[7] S. Granger and X. Pennec, “Multi-scale em-icp: A fast and robust approach for surface registration,” in Proceedings of European
Conference on Computer Vision (ECCV), pp. 418–432, 2002.
[8] B. Jian and B. C. Vemuri, “A robust algorithm for point set registration
using mixture of gaussians,” in Proceedings of International Conference
on Computer Vision (ICCV), vol. 2, pp. 1246–1251, 2005.
[9] E. Hasanbelliu, L. S. Giraldo, and J. C. Principe, “Information theoretic
shape matching,” Transactions on Pattern Analysis and MachineIntelligence (PAMI), vol. 36, no. 12, pp. 2436–2451, 2014.
[10] F. Wang, B. C. Vemuri, and A. Rangarajan, “Groupwise point
pattern registration using a novel cdf-based jensen-shannon divergence,” in Proceedings of Conference on Computer Vision and Pattern
Recognition (CVPR), vol. 1, pp. 1283–1288, 2006.
[11] J. C. Gower and G. B. Dijksterhuis, Procrustes problems (Oxford
Statistical Science Series Book), vol. 30. 2004.
[12] S. Gold, A. Rangarajan, C.-P. Lu, S. Pappu, and E. Mjolsness,“New algorithms for 2d and 3d point matching: Pose estimation andcorrespondence,” Pattern Recognition, vol. 31, no. 8, pp. 1019–1031,
1998.
[13] D. Chetverikov, D. Svirko, D. Stepanov, and P. Krsek, “The trimmed
iterative closest point algorithm,” in Proceedings of Conference on
Computer Vision and Pattern Recognition (CVPR), vol. 3, pp. 545–
548, 2002.
[14] J. Yang, H. Li, D. Campbell, and Y. Jia, “Go-icp: A globally optimal
solution to 3d icp point-set registration,” Transactions on PatternAnalysis and Machine Intelligence (PAMI), vol. 38, no. 11, pp. 2241–
2254, 2015.
[15] M. Greenspan and M. Yurick, “Approximate kd tree search for
efficient icp,” in Proceedings of International Conference on 3-D
Digital Imaging and Modeling, pp. 442–448, 2003.
[16] K.-L. Low, “Linear least-squares optimization for point-to-plane icp
surface registration,” Tech. Rep. 10, Chapel Hill, University of North
Carolina, 2004.
[17] 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.
[18] A. Segal, D. Haehnel, and S. Thrun, “Generalized-icp,” in Proceedingsof Robotics: Science and Systems, vol. 2, p. 435, 2009.
[19] Z. Zhou, J. Zheng, Y. Dai, Z. Zhou, and S. Chen, “Robust non-rigid
point set registration using student’s-t mixture model,” PloS one, vol. 9,
no. 3, p. e91381, 2014.
[20] Z. Min, J. Wang, and M. Q.-H. Meng, “Robust generalized point
cloud registration using hybrid mixture model,” in Proceedingsof International Conference on Robotics and Automation (ICRA),
pp. 4812–4818, 2018.
[21] S. Amari, Information geometry and its applications (Applied
Mathematical Sciences Book), vol. 194. 2016.
[22] Y. Qiao and N. Minematsu, “A study on invariance of f -divergenceand its application to speech recognition,” Transactions on Signal
Processing (TSP), vol. 58, no. 7, pp. 3884–3890, 2010.
[23] P. Biber and W. Straßer, “The normal distributions transform: A new
approach to laser scan matching,” in Proceedings of International
Conference on Intelligent Robots and Systems (IROS), vol. 3, pp. 2743–
2748, 2003.

[24] T. Stoyanov, M. Magnusson, and A. J. Lilienthal, “Point set registration
through minimization of the l2 distance between 3d-ndt models,” inProceedings of International Conference on Robotics and Automation(ICRA), pp. 5196–5201, 2012.
[25] B. W. Silverman, Density estimation for statistics and data analysis.
2018.
[26] H. Chui and A. Rangarajan, “A feature registration frameworkusing mixture models,” in Proceedings of Workshop on Mathematical
Methods in Biomedical Image Analysis, pp. 190–197, 2000.
[27] W. Qiao et al., “Asymptotics and optimal bandwidth for nonparametricestimation of density level sets,” Electronic Journal of Statistics, vol. 14,
no. 1, pp. 302–344, 2020.
[28] J. Jiang, Y.-L. He, D.-X. Dai, and J. Z. Huang, “A new kernel density
estimator based on the minimum entropy of data set,” Information
Sciences, vol. 491, pp. 223–231, 2019.
[29] H. Chui and A. Rangarajan, “A new point matching algorithm
for non-rigid registration,” Journal on Computer Vision and Image
Understanding (CVIU), vol. 89, no. 2-3, pp. 114–141, 2003.
[30] Z. J. Yew and G. H. Lee, “Rpm-net: Robust point matching using
learned features,” in Proceedings of Conference on Computer Vision
and Pattern Recognition (CVPR), 2020.
[31] D. Campbell and L. Petersson, “An adaptive data representation
for robust point-set registration and merging,” in Proceedings ofInternational Conference on Computer Vision (ICCV), pp. 4292–4300,
2015.
[32] L. Peng, G. Li, M. Xiao, and L. Xie, “Robust cpd algorithm for nonrigid point set registration based on structure information,” PloS one,
vol. 11, no. 2, 2016.
[33] J. Zhu, C. Jin, Z. Jiang, S. Xu, M. Xu, and S. Pang, “Robust point
cloud registration based on both hard and soft assignments,” Optics &Laser Technology, vol. 110, pp. 202–208, 2019.
[34] J. Ma, J. Zhao, J. Jiang, and H. Zhou, “Non-rigid point set registration
with robust transformation estimation under manifold regularization,”in Proceedings of AAAI Conference on Artificial Intelligence, 2017.
[35] R. Sinkhorn, “A relationship between arbitrary positive matrices anddoubly stochastic matrices,” Annals of Mathematical Statistics, vol. 35,
no. 2, pp. 876–879, 1964.
[36] B. Curless and M. Levoy, “A volumetric method for building complex
models from range images,” in Proceedings of SIGGRAPH, pp. 303–
312, 1996.
[37] S. Choi, Q.-Y. Zhou, S. Miller, and V. Koltun, “A large dataset ofobject scans,” arXiv preprint arXiv:1602.02481, 2016.
[38] F. Pomerleau, M. Liu, F. Colas, and R. Siegwart, “Challenging data
sets for point cloud registration algorithms,” International Journal of
Robotics Research (IJRR), vol. 31, no. 14, pp. 1705–1711, 2012.