1992 - A Method For Registration of 3-D Shapes - ICP - OK

IEEE TRANSACFIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 14, NO.
2, FEBRUARY 1992
239
A Method for Registration of 3-D Shapes

Paul J. Besl, Member, IEEE, and Neil D. McKay
AbstractThis paper describes a general-purpose, representation-independent method for the accurate and computationally efficient registration of 3-D shapes including free-form curves and surfaces. The method handles the full six degrees of freedom and is based on the iterative closest point (ICP) algorithm, which requires only a procedure to nd the closest point on a geometric entity to a given point. The ICP algorithm always converges monotonically to the nearest local minimum of a meansquare distance metric, and experience shows that the rate of convergence is rapid during the rst few iterations. Therefore, given an adequate set of initial rotations and translations for a particular class of objects with a certain level of shape complexity, one can globally minimize the mean-square distance metric over all six degrees of freedom by testing each initial registration. For example, a given model shape and a sensed
data shape that represents a major portion of the model shape
general, uni ed approach, which generalizes to n dimensions
and provides solutions to 1) the point-set matching problem

without correspondence and 2) the free-form curve matching problem. The algorithm requires no extracted features, no
curve or surface derivatives, and no preprocessing of 3-D data,

except for the removal of statistical outliers.
The main application of the proposed method as described

here is to register digitized data from un xtured rigid objects
can be registered in minutes by testing one initial translation and a relatively small set of rotations to allow for the given level of model complexity. One important application of this method is to register sensed data from un xtured rigid objects with an ideal geometric model prior to shape inspection. The described method is also useful for deciding fundamental issues such as the congruence (shape equivalence) of different geometric representations as well as for estimating the motion between point sets where the correspondences are not known. Experimental results show the capabilities of the registration algorithm on point sets, curves, and surfaces. Index Terms Free-form curve matching, free-form surface matching, motion estimation, pose estimation, quaternions, 3-D registration.
with an idealized geometric model prior to shape inspection. When inspecting shapes using high-accuracy noncontact measurement devices [4] over a shallow depth of eld, the uncertainty in different sensed points does not vary by much. Therefore, for purposes of simplicity and relevance to inspection applications based on thousands of digitized points, the case of unequal uncertainty among points is not considered. Similarly, the removal of statistical outliers is considered a
preprocessing step, is probably best implemented as such,
and will also not be addressed. In the context of inspection

applications, the assumption that a high-accuracy noncontact
measurement device does not generate bad data is reasonable since some sensors have the ability to reject highly uncertain measurements. The proposed shape registration algorithm can be used with the following representations of geometric data: 1) Point sets 2) line segment sets (polylines) 3) implicit curves: '(a:,y,z) = 0
4: parametric curves:
I. INTRODUCTION
GLOBAL AND local shape matching metrics for freeform curves and surfaces as well as point sets were described in [3] in an attempt to formalize and unify the description of a key problem in computer vision: Given 3D data in a sensor coordinate system, which describes a data
I triangle sets (faceted surfaces)

I implicit surfaces: g($,y,z) = 0
P} E! xax \
| parametric surfaces: (a:(u,v),y(u,"v),z(u,'v)).
This covers most applications that would utilize a method
shape that may correspond to a model shape, and given a

model shape in a model coordinate system in a different
to register 3-D shapes. Other representations are handled by providing a procedure for evaluating the closest point on the
given shape to a given digitized point.
geometric shape representation, estimate the optimal rotation

and translation that aligns, or registers, the model shape and
the data shape minimizing the distance between the shapes

and thereby allowing determination of the equivalence of the shapes via a mean-square distance metric. Of key interest to
This paper is structured as follows: Several relevant papers from the literature are rst reviewed. Next, the mathematical preliminaries of computing the closest point on a shape to a
given point are covered for the geometric representations men-
tioned above. Then, the iterative closest point (ICP) algorithm

is stated, and a theorem is proven concerning its monotonic convergence property. The issue of the initial registration states is addressed next. Finally, experimental results for point sets, curves, and surfaces are presented to demonstrate the
many applications is the following question: Does a segmented

region from a range image match a subset of B-spline surfaces
on a computer-aided-design (CAD) model? This paper provides a solution to this free-form surface matching problem
as de ned in [3] and [5] as a special case of a simple,

Manuscript received October 30, 1990; revised May 6, 1991. The authors are with the Computer Science Department, General Motors Research Laboratories, Warren, MI 48090-9055. IEEE Log Number 9102686.
capabilities of the ICP registration algorithm. II. LITERATURE REVIEW

Relatively little work has been published in the area of regis-
tration (pose estimation, alignment, motion estimation) of 3-D
0162-8828/92$03.00 1992 IEEE
240
IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 14, NO. 2, FEBRUARY 1992
free-form shapes. Most of the existing literature addressing global shape matching or registration have addressed limited
classes of shapes, namely, 1) polyhedral models, 2) piecewise-
(super)quadric models [2], [27], and 3) point sets with known

correspondence. The reader may consult [6] and [14] for pre-
1985 work in these areas. For a sampling of other more recent related work not discussed below, see [8], [10], [12], [13],
by investigating shape matching based on generalized shape polynomials. This demonstrated some interesting theoretical results but remains to be demonstrated for practical use on complex surfaces. Szeliski [54] also describes a method for estimating motion from sparse range data without correspondence between the points and without feature extraction. His primary goal was
to create a method for estimating the motion of the observer
[19], [20], [24], [26], [34], [35], [37]. [39], [44]. [46], [48],
[53], [58], [59].
Historically, free-form shape matching using 3-D data was
done earliest by Faugeras and his group at INRIA [18], where

they demonstrated effective matching with a Renault auto part
(steering knuckle) in the early 1980s. This work popularized the use of quaternions for least squares registration of corresponding 3-D point sets in the computer vision community. The alternative use of the singular value decomposition (SVD)
algorithm [23], [1], [49] was not as widely known in this time
between two range image frames of the same terrain. Given the set of points from one frame, he applies a smoothness assumption to create a smoothing spline approximation of the points. Then, a conventional steepest descent algorithm is used to rotate and translate the second data set so that it minimizes the sum of the covariance-weighted z differences between the points and the surface. His approach is based on a regular my-grid structure, and true 3-D point-to-surface distances are
not computed. The steepest-descent approach is a slower
frame. The primary limitation of this work was that it relied on the probable existence of reasonably large planar regions
within a free-form shape.
alternative to reaching the local minima than our proposed

ICP algorithm described below. Szeliski uses optimal Bayesian
Schwartz and Sharir [50] developed a solution to the freeform space curve matching problem without feature extraction
mathematics to allow him to downweight noisier values at longer ranges from a simulated range nder. For navigation
range imaging sensors, the uncertainty in data points vary
in late 1985. They used a nonquaternion approach to computing the least squares rotation matrix. The method works well
signi cantly from the foreground to the background. For highaccuracy sensors with shallow depths of eld, the uncertainty
with reasonable quality curve data but has dif culty with very noisy curves because the method uses arclength sampling of
the curves to obtain corresponding point sets.
variation between points is orders of magnitude less and is of much less concern. Szeliski provides experimental results
for synthetic terrain data and a block. The terrain data motion
Haralick er al. [28] addressed the 3-D point-set pose estimation problem using robust methods combined with the
test was a simple translation along one axis: a 1-D correlation

problem. His block test did involve six degrees of freedom,
least squares SVD registration approach, which provided a

robust statistical alternative to the least squares quaternion or SVD point set matching. This algorithm is able to handle
but the block is a very simple shape. Overall, this work

presents some interesting ideas, but the experimental results are unconvincing for applications.
statistical outliers and could theoretically be substituted for

our quatemion-based algorithm as long as the determinant of
Horn and Harris [33] also addressed the problem of estimating the exact rigid-body motion of the observer given sequentially digitized range image frames of the same terrain. They describe a range rate constraint equation and an elevation
the orthonormal matrix is strictly a positive one. A recent

conference proceedings [47] contains new contributions on this subject. Horn [31] derived an alternative formulation of Faugerass method [18] of least squares quaternion matching that uses
rate constraint equation. The result is a noniterative least squares method that provides a six-degree-of-freedom motion
estimate as long as the motion between frames of data is
the maximum eigenvalue of a 4 >< 4 matrix instead of the

minimum eigenvalue. Horn [30] and Brou [11] also developed
relatively small. This method is much quicker than the one

proposed by Szeliski, but it is not clear that this method generalizes to arbitrary rotations and translations of a shape.
the extended Gaussian image (EGI) methods allowing the

matching of convex and restricted sets of nonconvex shapes
Kamgar-Parsi er al. [36] also describe a method for the

registration of multiple overlapping range images without distinctive feature extraction. This method works very well using the level sets of 2.5-D range data but is essentially restricted to the three degrees of freedom in the plane since the
based on surface normal histograms. Taubin [55] has done some interesting work in the area of
implicit algebraic nonplanar 3-D curve and surface estimation
with applications to position estimation without feature extraction. He describes a method of approximating data points with implicit algebraic forms up to the tenth degree using
work was addressed toward piecing together terrain map data.

Li [38] addressed free-form surface matching with arbitrary
an approximate distance metric. Global shapes (not occluded

shapes) can be identi ed based on generalized eigenvalues, and
rotations and translations. His method forms an attributed

relational graph of fundamental surface regions for data and model shapes and then performs graph matching using an
the registration transformation can be recovered. The method

is shown to be useful for complete planar curve and space
inexact approach that allows for variability in attributes as

well as in graph adjacency relationships. This seems to be
curve shapes, but it is unclear that the effectiveness generalizes

well to more complicated surfaces, such as terrain data or a
a reasonable approach but relies on extraction of derivativebased quantities. Experimental results are shown for a coffee cup and the Renault auto part; see also Wong er al. [60] for other related work using attributed graphs for 3-D matching.
human face. Taubin has stated that the numerical methods of the approximate distance t tend to break down above the tenth
degree. He later [56] extended his work in shape description
BESL AND MCKAY: METHOD FOR REGISTRATION OF 3-D SHAPES
24]
The work of Gilbert and Foo [21] and Gilbert et al. [22] is related in that it addresses the computation of distance between two object shapes. Such methods could be the basis for similar shape matching techniques as are described below.
The major inconvenience with their method, though, is that
A. Point to Parametric Entity Distance In this section, a parametric curve and a parametric surface are treated as a single parametric entity 1F'(ri,'), where 11' = 'u. E R1 should be substituted for parametric curves, and ii : (u, '0) E R2 should be substituted for parametric surfaces
object shapes must be decomposed into convex subbodies, which is a problem that, in general, is not trivial for arbitrary
CAD models or for digitized 3-D data.
("R denotes the real line). The evaluation domain for a curve
is an interval, but the evaluation domain for a surface can be
any arbitrary closed-connected region in the plane. For more

information on parametric entities, such as Bezier and B-spline
III. MATHEMATICAL PRELIMINARIES

In this section, methods for computing the closest point to
curves and surfaces, see [9], [15][17], [42], [52].

The distance from a given point 13 to a parametric entity E is
a given point on the various geometric representations listed

above are described. First, the basic geometric entities are covered, followed by parametric entities, and, nally, implicit
dtri E) = r("u.)EE -1131 4051 F'(17))-
(6)
entities. The reader might consult Mortenson [42] on some of the items below for additional information. The Euclidean distance d('.F'1,'F2) between the two points F1 = (56112/1,21,) and F2 = ($2.3/2,Z2) is d('F1J-$2) I HT1 Tall I \/(372 11)2 + (3/2 ?J1)2 l' (Z2 - Z1)2- Lei A be a point set with Na points denoted 5,: A = {G5,} for 2' = 1, . . . , Na. The distance between the point 15 and the point set A is
The computations to compute the distance are not closed form

and are relatively involved. One method for computing point-
to-curve and point-to-surface distances is described below.

Sets of parametric entities are again straightforward once the
distance metric for an individual entity is implemented. Let

F be the set of Ne parametric entities denoted E,, and let
F=
for 2' : 1, Ne. The distance between a point 13' and
the parametric entity set F is
d(r?l A) = i{jn_ir1Na} d(1Ifi.-)-
(1)
d(r, F) ,{1j_1}N} d(r1, Er)
(7)
The closest point ti, of A satis es the equality d(g3, cifj) : d(g5', A). Let l be the line segment connecting the two points F1 and
F2. The distance between the point F and the line segment Z is
I H1111
"u;-I-11:1
The closest point 3], on the parametric entity set F satis es the equality d(p', 3],-) = d(p,
Our rst step towards computing the distance from a point to
a parametric entity is creating a simplex-based approximation

(line segments or triangles). For a parametric space curve
C : {F'('u.)}, one can compute a polyline L(C, 6) such that the

+ UF2
piecewise-linear approximation never deviates from the space
where u [0,1] and "v E [0,1]. The required closed-form computations are straightforward. Let L be the set of N; line
curve by more than a prespeci ed distance 6. By tagging each point of the polyline with its corresponding u argument values
of the parametric curve, one can obtain an estimate of the 'u,,,
segments denoted l,, and let L :
for z' : 1,. . . , N1. The
argument value of the closest point from the line segment

set.
distance between the point 15' and the line segment set L is
do: L) = {;13t;3N1}d<a.1a>~
equality d(p', ',-) d(13, L).
<3)
Similarly, for a parametric surface S : {'F'(u,'u)}, one can compute a triangle set T(S,6) such that the piecewisetriangular approximation never deviates from the surface by
The closest point 3}, on the line segment set L satis es the Let t be the triangle de ned by the three points F1 =
(SL'1,y1,/2'1), 772 I (.5l3g,_t/2,22), and F3 I (Ll5'3,f_lj3,Z3). Tll
more than a prespeci ed distance 6. By tagging each triangle

vertex with the corresponding (uyu) argument values of the parametric surface, one can obtain an estimate (ua, pa) of the argument values of the closest point from the triangle set. As a
distance between the point 13 and the triangle t is

d(;rT, t) = u+rgi3:1 ||uF1 + 111*; + wf}, 35']| (4)
result of these curve and surface procedures, one can assume that an initial value 1].], is available such that r"('&.',,) is very
close to the closest point on the parametric entity. The point-to-parametric-entity distance problem is ideal for employing a pure Newtons minimization approach when a
where u E [0, 1], o G [0,1], and w E [0,1]. The required
closed-form computations are again straightforward. Let T be the set of N, triangles denoted t,~, and let T : {t,-} for 21 = 1,..._,N,. The distance between the point 15' and the
triangle set T is given by
reliable starting point ii, is available. The scalar objective

function to be minimized is
f( Qt ) = ||t"(I1')r5'l|2Let V I [3/6911'] I be the vector differential
(8)
gradi-
4(5) T) = ,{j11_iI}Nt} d(r1ir)4051 17;") = d(r51 T)-
(5)
The closest point y on the triangle set T satis es the equality
ent operator (where t implies vector transpose). The minimum of f occurs when V f = 0. When the parametric entity is a surface, the 2-D gradient vector is V f = [f,,, f,,]*, and the 2-D Hessian matrix
242
is
quadratic
objective
function
subject
to
nonlinear
We el ii: iii]
by Mil) f.)(i-1) fm.(t7) f....(ii) fUU(,J)
f,,,,. The
constraint Minimize f('F') = ||'.F'p']|2 where '('F') : 0.

Lagrange multiplier system of equations [40]:
(18)
where the partial derivatives of the objective function are given
One approach to this problem is to form the augmented
= 2'Ff.(11') ('F(11') ~ 25) = 2ti(t1')('F'('J) 15') = 2t"i...(t1)(t*(t1') - 15') + 2Fi,(11')Fu('J) = 21*..(11')(t( ti) + 2Fi('@1)t%(iTT) = 21*i...(t1')(F( Eifi I5) + 2rZ(1?)Ft(t1')update formula for either
(10) vne+Yvao-0 an (11) an =0 (12) (13) where V 6/6f]* and solve this system of nonlinear equations via numerical methods. The number of equations (14) and unknowns for the nonlinear system is three for planar
curves, four for surfaces, and ve for implicitly de ned space entity
The curve case requires only computation of f,, and

Newtons
curves. Continuation methods [41] can be used to solve this

problem for algebraic entities even without a good starting
is
tnI=ar-Wv%naar*vnnJ
an
point, but a good starting point will allow the use of faster
methods, such as the multidimensional Newtons root nding
where rig = 11,. When using the starting point selection
method described above based on a simplex approximation

with a reasonably small 6, Newtons method for computing
method. From a numerical point of view, the parametric methods are much easier to deal with. From an applied
point of view, no industrial CAD systems store free-form
the closest point generally converges in one to ve iterations

and typically in three. The computational cost of Newtons
curves or surfaces in implicit form. For this reason, implicit

surfaces of interest are dealt with in our implemented system
method is very low in contrast with points.

B. Point to Implicit Entity Distance
nding good starting
either via special case mathematics (e.g., spheres) or via

a parametric form. Of course, if there were an application
where it was necessary to handle free-form implicit entities

in their implicit form [51], the above algorithm could be implemented. Taubin [55] uses an approximate distance algorithm that implies a simple update formula for surfaces and planar curves when g(F0) is nearly zero:
An implicit geometric entity is de ned as the zero set of
a possibly vector-valued multivariate function (F') : 0. The

distance from a given point 15' to an implicit entity I is
(K2521) = s(1")=0 311111 d(15'F) = 9("")=0 P1111 |lI" rill
(16)
_.
The calculations to compute this distance are also not closed form and are relatively involved. One method for computing point-to-curve and point-to-surface distances is outlined below.
I it ' m"
5{3L=5-11:;
_.
V9( it
(Ft)
<2)
This method is only exact if the in nite line with the direction
Vg(F) at the starting point "F0 intersects the implicit entity at a
Sets of implicit entities are straightforward once the distance

metric for an individual entity is implemented. Let J be the
point where the normal vector has that same direction. This is
not true in general, and the approximation is generally worse the further the point is from the implicit entity. Therefore,
set of NI parametric entities denoted I), and J = {Ik} for

is = l, N1. The distance between a point 15' and the implicit entity set J is
this result cannot be used if precise distance results are

required. C. Corresponding Point Set Registration
aao:%g@Mgaa>
an
All closest point (minimum distance) algorithms have been

The closest point 37,- on the implicit entity I, satis es the equality d(15',',-) = d(ji',.]). mentioned in forms that generalize to n dimensions. One more
necessary procedure for yielding the least squares rotation and

translation is reviewed. For our purposes, the quaternion-based
Our rst step towards computing the distance from a point to an implicit entity is creating a simplex-based approximation
(line segments or triangles) as was done for parametric entities
algorithm is preferred over the singular value decomposition

(SVD) method in two and three dimensions since re ections
[7]. Computing the point-to-line set or point-to-triangle set

distance yields an approximate closest point Fa, which can be used to compute the exact distance. The implicit entity distance problem is quite different from
are not desired. The SVD approach, based on the crosscovariance matrix of two point distributions, does, however,
generalize easily to n dimensions and would be our method of

choice for rt > 3 in any rt-dimensional applications. The basic solution of Horn [31] is described below, although the method
the parametric entity case where unconstrained optimization suf ces. To nd the closest point on an implicit entity
de ned by '('F') = 0 to a given point p, one must solve a constrained optimization problem to minimize
of Faugeras [18] is equivalent. Our summary stresses the role

of the SVD cross-covariance matrix, which is an important relationship not discussed in other work.
243
The unit quaternion is a four vector (IR = [q0q1q2q3]*, where

(10 2 0, and qg + qf + qg + qg = 1. The 3 x 3 rotation matrix
IV. THE ITERATIVE CLOSEST POINT ALGORITHM Now that the methods for computing the closest point on a geometric shape to a given point and for computing a
least squares registration vector have been outlined, the ICP
generated by a unit rotation quaternion is found at the bottom of this page. Let 51- = [q4q5q6] be a translation vector. The complete registration state vector rj is denoted rj [rj'R|rjT]. Let P = be a measured data point set to be aligned with
a model point set X = {:i:',}, where N,, = Np and where each point 5, corresponds to the point 5, with the same index. The
algorithm can be described in terms of an abstract geometric

shape X whose internal representation must be known to
execute the algorithm but is not of concern for this discussion.

Thus, all that follows applies equally well to 1) sets of points, 2) sets of line segments, 3) sets of parametric curves, 4) sets
mean square objective function to be minimized is
of implicit curves, 5) sets of triangles, 6) sets of parametric
fa) =
Z Ila - Rena - a-||*\

1 1
NP
<22)
surfaces, and 7) sets of implicit surfaces.
In the description of the algorithm, a data shape P is

moved (registered, positioned) to be in best alignment with
The center of mass /1}, of the measured point set P and the center of mass x for the X point set are given by
a model shape X. The data and the model shape may be

represented in any of the allowable forms. For our purposes,
the data shape must be decomposed into a point set if it is

not already in point set form. Fortunately, this is easy; the
Hp =
'_B2
IIM3
$511 and /I5-3 =
I-182
.;MF
S-. S
(23)
points to be used for triangle and line sets are the vertices
and the endpoints, and if the data shape comes in a surface or
The cross-covariance matrix Em of the sets P and X is given
curve form, then the vertices and endpoints of the triangle/line approximation (as described above) are used. The number
of points in the data shape will be denoted NP. Let N,, be
W
310$ =
'_U2
the number of points, line segments, or triangles involved in
(24) The cyclic components of the anti-symmetric matrix A,_,- = (Em ZIg:,,),,- are used to form the column vector A =
[A23 A31 A1g]T. This vector is then used to form the symmetric 4 X 4 matrix Q(E,,,,)
[(1% /1-p)(@11: /1-alt] = 3; ZZ;lPt$il Hatti.3 .M5
"0
-5
the model shape. As described above, the curve and surface closest-point evaluators implemented in our system require a framework of lines or triangles to yield the initial parameter values for the Newtons iteration; therefore, the number N, is
still relevant for these smooth entities but varies according to

the accuracy of the approximation. -0 The distance metric d between an individual data point p and a model shape X will be denoted
d(i5ZX) =-" :rEX 11i11||i"15'||where I3 is the 3 >< 3 identity matrix. The unit eigenvector
q}} = [qg q1 Q2 q3]* corresponding to the maximum
(23)
The closest point in X that yields the minimum distance is

denoted 3] such that d(j5',g,T) : d(p',X), where '57 E X. Note
eigenvalue of the matrix Q(E,,,) is selected as the optimal

rotation. The optimal translation vector is given by
that computing the closest point is O(N,,) worst case with expected cost log When the closest point computation
(from 13 to X) is performed for each point in P, that process is
worst case O(NpN,,). Let Y denote the resulting set of closest
(IT = tie R(TR)/In:
(26)
points, and let C be the closest point operator:
This least squares quaternion operation is O(N,,) and is denoted as
Y=qRxy
am
(til dms) = Q(P, X)
(27)
Given the resultant corresponding point set Y, the least squares registration is computed as described above:
where rim, is the mean square point matching error. The nota-
(ii <1) = Q(P= Y)
(39)
tion rj'(P) is used to denote the point set P after transformation

by the registration vector Q.
The positions of the data shape point set are then updated via P = rj'(P).
2 2 _ <10 '1' Q1 " (1% (Iii
22(q1q2 - qoqal
2(q1q3 + (I092)
R=
2@m+mm)
2(q1q3 qoq2)
%+
2(q2qa + qoqil
2@m~mm)
<13 + qg qj qg
(3)
244
A. ICP Algorithm Statement

The ICP algorithm can now be stated:
' The point set P with Np points from the data shape
would yield a smaller mean square error than the least squares registration, which cannot possibly be the case. Next, let the
least squares registration rj}, be applied to the point set P0, yielding the point set P;,+1. If the previous correspondence to
and the model shape X (with N, supporting geometric primitives: points, lines, or triangles) are given.
' The iteration is initialized by setting P0 = P, (I0 =
the set of points Yk were maintained, then the mean square

error is still dk, that is
[1,0,0,0,0,0,0]* and k = 0. The registration vectors

are de ned relative to the initial data set P0 so that the nal registration represents the complete transformation. Steps 1, 2, 3, and 4 are applied until convergence within a tolerance 'r. The computational cost of each operation is given in brackets. a. b. Compute the closest points: Yk : C (Pk, X) (cost:
dk =
i_i G2
$3.2 We -i5,k+1||2-
(33)
However, during the application of the subsequent closest

point operator, a new point set Y;,,+1 is obtained: Y;,+1 I C(P;,+1,X). It is clear that
0(NpN,,) worst case, 0(Np log N,,,) average).

Compute the registration: (fIl==d;,) = Q(P0,Y;,)
|l?Ii.i=+1 - Z7t,k+1|| S lltjtk -}7t,k+1l| fol @3611 i = 11NP (34) because the point 311;, was the closest point prior to transformation by rj}, and resides at some new distance relative to 13',;,;,+1.
If g]',;,;,+1 were further from g'5',,;,+1 than Qk, then this would
c.
d.
(cost: O(Np)). Apply the registration: P;,,+1 : q";',(P0) (cost:
0(Np))Terminate the iteration when the change in meansquare error falls below a preset threshold *r > 0 specifying the desired precision of the registration:
dk dk+1 < '7'.
directly contradict the basic operation of the C closest point operator. Therefore, the mean square errors ek and dk must obey the following inequality:
0 5 d,.+1 g e;,+1 g at g 6,, for all k.
(35)
If a dimensionless threshold is desired, one can replace r

with rt/tr(Z,,), where the square root of the trace of the covariance of the model shape indicates the rough size of the model shape. B. Convergence Theorem A convergence theorem for the ICP algorithm is now stated
The lower bound occurs, of course, since mean-square errors

cannot be negative. Because the mean-square error sequence is nonincreasing and bounded below, the algorithm as stated
above must converge monotonically to a minimum value. Q.E.D. Experimentally, we nd fast convergence during the rst
few iterations that slows down as it approaches the local
minimum. Even at this slow pace, somewhere between 30

and S0 iterations yields excellent results: dk 2: 0.1% of model
and proved. The key ideas are that 1) least squares registration generically reduces the average distance between corresponding points during each iteration, whereas 2) the closest point
shape size. The convergence can be accelerated using a simple

additional operation described in the next section. C. An Accelerated ICP Algorithm
determination generically reduces the distance for each point

individually. Of course, this individual distance reduction also
reduces the average distance because the average of a set of smaller positive numbers is smaller. We offer a more elaborate
explanation in the proof below. Theorem: The iterative closest point algorithm always converges monotonically to a local minimum with respect to the
The accelerated ICP algorithm uses a minor variation on the basic line search methods of multivariate unconstrained
minimization [45]. As the iterative closest point algorithm
proceeds, a sequence of registration vectors is generated: (I1,

(jg, (Y3, 6'4, (I5, (I6, . . ., which traces out a path in the registration
mean-square distance objective function.

Proof. Given Pk = {13',;,,} = rj},(P0) and X, compute the
state space from the identity transformation toward a locally

optimal shape match. Consider the difference vector sequence de ned by
set of closest points Yk :
as prescribed above given
the internal geometric representation of X. The mean squared error e;, of that correspondence is given by
Alli I (Iii "- i~1
(36)
1 N _ _ Bk I *" Z!/tr *Ptk||2Np 1 1
which de nes a direction in the registration state space. Let (31)

the angle in 7 space between the two last directions be denoted
91 I CO5 1 Em
Afli<tA(Il-1
The Q operator is applied to get Q}, and dk from that corre-
spondence:
.N
||Aqt||llAqt_1l|
(37)
and let 66 be a suf ciently small angular tolerance (e.g., 10).
If
P % "-F D
dk = -A7 Z llytt R(qkR)PvI0 qkT||2~

p 1 1
(32)
Elk < and l9;,,_1 < (38)
It is always the case that dk ek. Suppose that dk > ek. If this were so, then the identity transformation on the point set
then there is good direction alignment for the last three registration state vectors: Q,-1,, }.,_1, and rj},_2. Let dk, dk_1,
BESL AND MCKAY: METHOD FOR REGISTRATION OF 3-D SHAPES 245

__ _ *'_'*' ' _"_'H"_ z; LIT
Linear Approximation Parabola
The Accelerated I PAl 'h C gem m
J
I .
-<s>
m=".::r....*"..=:....
\.
Mean Square Em Axis
d(i-2) d-1 ~._\
\____1~mmumBmr
q(i-2 in
Q
qti-3)
iv
q(1)
8e(1_I)
Linear Update
I l
'
:
Plane Represents 7-D Space
_ .
-Ems w _ ,_,,,_A_,, ___r __ ._
\ Tz
A-_ ._ _
A ,
J
mo
: Tx
, _ ,
g
_ 1' _-_l~-1-ig>4_.
,_
58(1)
Parabola Update
Z0-1'-I>l"0lZ>W"|
AL.///l
Fig. 1.
Consistent direction allows acceleration of the ICP algorithm.
and d;,_2 be the associated mean square errors, and let '01,,
v;,_1, and v;.,_;_> be associated approximate arc length argument

values:
B 1
ll)
l ; L . 5
'5
r ; *1
,::,
eido
f
I l
ql
'---_
I Jm __ ___ _ 1 ll 1_ l
zo-
l __ _ _"C i_ if i__'I - C _ Cjlm T 7 *-
4-4- _.-i --'7
CumuIltivcArcLtglh
" - 1 -%==' Eual -
Uh = U vt-1 = -|lAli<||,va-2 = "l|A(lit-1|] + 'vxt~1- (39) See Fig. 1 for a picture of the situation. Next, a linear
approximation and a parabolic interpolant to the last three
Fig. 2. Various quantities plotted against iteration count for the basic nonaccelerated ICP algorithm.
data points are computed:

d1(o) = alt; + bl d2(o) = 0.2212 + b2U + C2 (40)
which gives us a possible linear update, based on the zero crossing of the line, and a possible parabola update, based on
the extremum point of the parabola:
U1 = -01/(11 > 0 U2 I -02/20,2.
using the new registration with the last two steps and move to the appropriate minimum. This has not been necessary in our experience. To be rigorous, one can simply ignore the suggested update if it causes a worse mean square error. To give a quantitative example comparison, the registration
values, RMS error, maximum error, angular change, and cumulative arc length values were recorded during 50 iterations
To be on the safe side, we adopt a maximum allowable value omax. The following logic is used to perform an attempted update: 1)If0<Ug <'u1 <'vm,,Xor0 < '02 < vmax < "01, use the parabola-based updated registration vector: (fl, = (Ir + 'ugAcj},/|[Atj';,|| instead of the usual vector if}, when performing the update on the point set, i.e.,
of both the basic and accelerated ICP algorithms during the same free-form surface matching test. The results for the basic ICP algorithm are shown in Fig. 2. Note the smooth character
of all the curves. The most important feature is that the cos(6t9)
plot indicates a consistent direction of updates for all but the

rst few iterations. In contrast, the accelerated ICP algorithm shows the desirable jumpy behavior as seen in Fig. 3. In
Pa+1 = <I'i,(P0)-
addition, note how most quantities get close to their nal values after the rst acceleration step and very close after two. The acceleration steps occur whenever a V-shaped dip occurs
in the plot of cos(6l9) versus the iteration count. D. Alternative Minimization Approaches The ICP algorithm allows us to move from a given starting
2) If0<v1 <v2<vm.-ix 0r0<"v1 <vI11&X <v20r

U2 < 0 and 0 < o1 < omax, use the line-based updated
registration vector c = (Tr + o1Aq'},/HArj';,|| instead of
the usual vector if}.

3) If both v1 > omax and '02 > vmax, use the maximum allowable update (T2, = iii + t1m,,,,,Atj;,/||Arj';,|| instead of the usual vector cj},,. We have found experimentally that setting vmax : 25|]Atj'k|] adaptively has provided a good sanity check on the updates allowing the iterative closest point algorithm to move to the
point to a local minima in 7 space relatively quickly in

comparison with other possible alternatives. Each iteration
requires only one evaluation of the closest point operator:

the most expensive computation. Any optimization method
local minimum with a given degree of precision in many fewer

steps. A nominal run of more than 50 basic ICP iterations for a given value of 1' is typically accelerated to 15 or 20 iterations. If the updated registration vector were somehow to overshoot the minimum enough to yield a worse mean square error, it would be advantageous to construct a new parabola
that does not use explicit vector gradient estimates, such as Powells direction set method, the Nelder-Mead downhill
simplex method, or simulated annealling, requires literally hundreds to tens of thousands of closest point evaluations. These numbers are based on tests done to simulate the action of the least squares registration step involved in one ICP iteration but using instead Powells direction set method and the Nelder-Mead method from [45].
246
I L
/\r v
V
ytos (69)
v V 1:
Horizontal Axis: Number of Iterations
has a bad effect on most nonlinear optimization algorithms.
To summarize, any method that allows one to move from

an initial state to its corresponding local minimum could
XlI
Um Error
theoretically be used in place of the ICP algorithm. For

ES Error l. _. .
example, consider Szeliskis [54] work with steepest descent and three rotation angles. However, the arguments above indicate that one would have a hard time trying to nd an
algorithm that was only ten times slower on average. The key
bene t of the ICP algorithm is that the convergence is fast and monotonic. No expensive closest point evaluations are spent
><*< N
on registration vectors that have worse mean square errors than

Z
the current state. Because of the ICP convergence theorem, one
does not have to feel around in the multidimensional space

to determine the direction in which to move.
V.
2Q-"Cl>)-%]C)>['F<gU;Zj>7J'-l in
THE SET OF INITIAL REGISTRATIONS
> -
Cumulative Arc Length

\% Ex & ,9) -2 hi -1
Even though the ICP algorithm must converge monotonically to a local minimum from any given rotation and translation of the data point set, it may or may not converge on the desired global minimum. The only way to be sure is to
Fig. 3. Various quantities plotted against iteration count for the accelerated ICP algorithm.
nd the minimum of all the local minima. The problem with

reaching the desired global minimum with certainty is that it
is dif cult to precisely characterize in general the partitioning

of the registration state space into local minima wells (regions
Any optimization method that uses explicit vector gradients, such as steepest descent, conjugate gradient, and variable
metric schemes, will require at least seven closest point
of attraction) because this partitioning is potentially different for every possible different shape encountered.
To be precise, consider a 6-D state space Q, where the
evaluations for each numerical gradient evaluation. Therefore,

such a method would have to converge in three or four
quaternion component qo is determined from the other quaternion components: qo = (/1 (qf + q + (11% The actual
iterations to be competitive with the accelerated ICP method.

Such generic methods generally require many more than three iterations with the number of required closest point evaluations running well over 100 even in ideal circumstances. If a pure
state space Q is a subset of the space Q = [-1, 1]? >< R3, where R : (oe,+oo) is the real line. The subset Q is
speci ed by the inside or on the unit 3 sphere constraint
that qf + qg + qg 5 1. Therefore, Q may be viewed as a type

of hyper-cylinder in 6 space.
numerical Hessian-based Newtons method were used, the numerical gradient and Hessian computations would require at least 13 closest point evaluations per iteration, implying
that the iteration would have to converge in two iterations
For any given nonpathological shape X that represents

a real-world surface or object (e.g., pathological shape de-
scriptions based on sin(1/:12) near zero not allowed) and for

any given point set Preg already correctly registered with
to be competitive with the accelerated ICP algorithm. A pure Newtons method might require only three iterations if the initial point were already well into the region of attraction surrounding a local minimum, but the initial iterations would
not be handled well by Newtons method. Whenever an accelerated parabolic update takes place after three basic ICP steps, we can get nearly quadratic convergence
X, consider that any initial state Q G Q of the point set P = q(P;@g) will converge to a local minimum as it is
matched to X. There are a nite number of local minima
Nm(X, P) after one has xed X and P. (The shape X is

considered pathological if this is not true.) Let \I1(X , P) be
the set of all local minima:
for less than steepest descent cost. This is an interesting

accomplishment for a function where derivatives cannot be
\I1tX. P) = {..}.t:1
This allows us to state that
<42)
evaluated. Note that the steepest descent gradient direction

is not deliberately computed; we merely observe when a
This induces a natural partitioning of Q into equivalence
classes, labeled \I1,,, where every value of if that converges

via the ICP algorithm to ipn is a member of the class \I1n.
Nm
consistent direction is being followed. Other problems involved with using general-purpose optimization methods are the following: 1) If any angles are
used as in [54], angular cycles across 360 must be handled correctly, and 2) if a unit quaternion becomes a nonunit quaternion, as would be expected taking arbitrary direction steps
ti: U qr, and \I1.,,r|\I1,.,, =irnm.

'n.=1
(43)
in 4 space, the quaternion must be renormalized somewhere.

Unfortunately, if the objective function evaluator changes the
Let 111 be the equivalence class that maps to the correct global minimum 1111. To guarantee that the global minimum is found for a given
values in the state vector during the optimization iteration, this
shape X and a given set P not already registered with X, one
247
must use an appropriate set of initial states so that transforming
2) transform it rst so that the centers of mass of P and X
P by at least one initial registration will place the point set into the correct equivalence class \I11 of registrations. This
allows it to converge to the correct global minimum 1/)1. Two
coincide, and then apply ICP. We have found no differences in the nal registration results between 1) not translating and
2) pretranslating the point set when using an adequately large
fundamental questions are 1) how to construct an initial set

of states for any given object that guarantees a correct global
set of initial rotations (e.g., 24). In fact, any translation state

suf ces because the ICP algorithm is very insensitive to the
minimum and 2) how to construct an initial set of states that

guarantees all shapes in a given class of shapes when those
initial translation state when used for global shape matching.

We have observed that pretranslating the data point sets center
shapes converge to the correct respective global minimum. By using a suf ciently dense uniform sampling of quaternions on the unit sphere combined with a suf ciently dense sampling of translation vectors occupying the total volume about the shape X, it is possible to determine the complete nite set of local minima with a suf ciently small probability of error for one given object. One could construct a set of
of mass to the model shapes center of mass generally saves a few iterations (e.g., 2 to 4) out of the usual 20 or so total. A further simpli cation in the global shape matching algorithm can be accomplished for most generic shapes, where principal moments demonstrate some level of distinctness. Let
pr Z pg Z pz be the square roots of the eigenvalues of Zlp, and let 11,, Z ry Z rz be the square roots of the eigenvalues
initial states to include all these local minima solutions along

with the halfway states between the h nearest neighbors (e.g.,
of Em. If the following sets of conditions hold: Py S 012%

ry agrm
h = 12) to attempt to avoid having a set of initial states that

lie on or near the boundaries between the equivalence classes.
Pi. S 042%
T2 O!2'l"y
(45)
(46)
A method that could be useful for computing guaranteed

initial states for a set of model shapes is to use a 6-D
occupancy array to compute hypervoxel-based descriptions

of the equivalence classes for each shape of interest and
for a speci ed value of 0:2, e.g., O52 = 1/\/2 z 0.71, one
then computing an overall partitioning by re ning the
rst
can reliably match the basic shape structure of data and model using only the eigenvectors of the matrices E1, and Em. Again,
the exact value of 0:2 could be computed for any given set of objects and any given level of sensor noise via exhaustive testing if needed. In this case of eigenvalue distinctness,
shapes partition by intersecting the equivalence classes of
each subsequent shape. Such methods can be very memory intensive; a 20 >< 20 >< 20 >< 20 >< 20 >< 20 hypercubic-hypervoxel
grid of the smallest hypercylinder containing all relevant
the identity transformation and the 180 rotations about the

eigenvector axes corresponding to 1",, ry, and T2 provide a very good set of only four initial rotations that will yield the correct global minimum for a wide class of model shapes. If two of the three eigenvalues are approximately equal but signi cantly different from the third for both data and model
registrations of all shapes of interest requires an 8-Mbyte

(64 Mb) array for a single object. Clearly, this is a test of
patience for anyone wishing to construct an initial set of states

customized to a given set of objects to yield guaranteed results,
but indeed, it is possible and relatively straightforward, except for the high dimension of the problem.
A. Initial States for Global Matching
shapes, the number of initial states need only be expanded for

rotations about the nonambiguous axis, thereby reducing the total number of initial rotational states.
If neither of the above cases for global matching hold true

(i.e., pa, z pg m pz and 1",, m ry re rz), then one must use a
There are simpler methods of dealing with the initial state

problem that are very effective on most shapes one comes across. Let us adopt the following de nition of the rst two
ne sampling of quaternion states that cover the entire surface

of the northern hemisphere of the unit 4 sphere uniformly. The rotation groups of the regular polyhedra, which have
moments of the distribution of geometry in P and X: /1}, =

P], and Z3, 2 E[(F ;1,,,)(F )ti',_,)t[F E X], where represents the sample expectation (averaging) operator. If the point data set P covers a signi cant portion of the model
El l e P]. /It = Elflf Xl En = E[(r5'/It-)(15'/Ip)|r5
been well known to crystallographers since the 1800s [29],

provide a convenient set of uniformly sampled initial rotations:
(a) 12 tetrahedral group states, (b) 24 octahedral/hexahedral

group states, and (c) 60 icosahedral/dodecahedral group states.
shape X such that the condition
The tetrahedral states are a proper subset (subgroup) of the
a1\/tr(E,,) 5 ,/it~(Ep) 5 \/tr(E,,)
(44)
holds for a suf ciently large factor, say 0:1 = 1/\/2 m 0.71, then we have found that it is generally not necessary to use
multiple initial translation states, as long as enough rotation

states are used. This factor (131 is the allowable occlusion percentage for global matching. The exact value of O51 could be computed for any given class of object shapes via exhaustive
octahedral/hexahedral states and the icosahedral/dodecahedral states. The octahedral/hexahedral states are not properly contained in the icosahedral/dodecahedral states. For a convenient listing of these rotations in quaternion form, see Appendix A of [32]. From an implementation point of view, one has the option
of using precomputed lists or nested loops. For the nested

loop case, all normalized combinations of qo = [1, 0}, ql = {+1,0, -1}, Q2 = {+1, 0, -1], and Q3 ={-l-1,0, -1} provide an easy-to-compute set of 40 rotation states that includes the
testing if that is desired.

There are two reasonable options for the initial translation
tetrahedral and the octahedral/hexahedral groups. (One must

ensure that the rst nonzero quaternion component is positive
state: 1) Apply the ICP algorithm directly to the point set P

using multiple rotation states about its center of mass /11,, or
to avoid duplication of states.) For a really complicated set
248
of shapes, all normalized combinations of Q9 =- {1,0.5,0},
The ICP algorithm is still useful for local matching problems

where the entire set of data points P matches a subset of
q, = {+1,0.5,0, -0.5, -1}, qg : {+1,0.5,0, -0.5, -1},

and (13 = {+1,0.5, 0, -0.5 l provides another easy-tocompute set of 312 initial rotation states. Another scheme for
the model shape X. The drawback is that more than one

initial translation must be used, which increases the cost of
a very dense sampling of states is to re ne the 60 states of

the icosahedral group by subdividing triangles using a 1 to 4 split and using an increased number of rotations about each axis speci ed by each vertex of the re ned icosahedron.
computing the correct registration. If N, is the number of

initial translation states, Nq is the number of initial rotation
states, and each closest point evaluation of P relative to X

costs O(N,,N$), in the worst case, the total cost of local matching is O(N,NqN,,N,,) as opposed to the cost of global matching O(NqN,,N,,). The number of initial translations is also dependent on the relative size of the data point set P
The general rule of thumb is the more complicated the

object, the more initial states required. Any method of getting a suf ciently dense, uniform distribution of quaternions over the northern hemisphere of the unit 4 sphere (or over the
compared with the model shape X. By de ning a quantity n

as the ratio of the sizes of the data and model shapes
full interior and surface of the unit 3 sphere) is adequate. In general, every application may want to use a customized set
of initial quaternions that will maximize the probability of
choosing a good starting point early for the shapes of interest.

B. A Counter Example
4=
21
(41)
where m(-) is a general measure that measures approximate

1) arc length if X is a curve, 2) area if X is a surface, and 3) volume if X is a volumetric point con guration, then one can estimate the basic qualitative behavior of the required number of translation states as N, m c(X)r7 for most shapes, where
Although the above methods for global shape matching will

work very well for many shapes with a small probability of
error, we can also state categorically that for any given xed
set of initial rotation states, one can construct a shape X that cannot be correctly registered by the algorithm. Begin with
c(X) is an approximate proportionality factor that depends on

the complexity of the shape X being matched.
Computing such an estimate m(X) is straightforward for

the shape X, but evaluating m(P) is more dif cult because P may simply be a point set, and the corresponding length, area,
a sphere of radius R. Add a thin spike of length S to the surface of the sphere for each speci ed axis and for each rotation about that axis as indicated by the fixed set of initial
rotation states. Next, add one or more spikes of length S + e anywhere that the extra spikes will t. If the extras will not t, make the original spikes thinner. Then, sample points on
and volume on the shape model X is unknown. One can use

a convex hull, surface Delaunay, or closest point connection
algorithm to get accurate measures for volumes, areas, and

arc lengths, but it is more likely that one would design an
the surface of this shape with any desired scheme so that there
is at least one point per spike, including one at the tip of each spike. The ICP algorithm combined with the given xed set
algorithm to tolerate up to a particular percentage of occlusion,

and a xed set of translations would be computed to handle
all objects in a given class of objects with up to that level of

occlusion. Examples of local shape matching are demonstrated in the next section.
VI. EXPERIMENTAL RESULTS
of initial rotation states will not be able to register a generic

repositioning of this point set with the original object in such a way that the longer (or shorter) spikes are correctly registered
with each other. It can safely be predicted that the proposed

registration algorithm will have dif culty correctly registering
sea urchins and planets. These shapes are characterized

as having almost exactly equal eigenvalues of the covariance matrix E, and as having small shape features at a ne scale relative to the overall shape. Of course, for any given xed set of object shapes, the set of initial rotations can be increased
This section is divided into three sections: 1) point set

matching, 2) curve matching, and 3) surface matching. All
programs were written in C. Any quoted approximate times are given for execution on a single-processor computer rated at 1.6
M ops on the 100 >< 100 double-precision Linpack benchmark. A. Point Set Matching In this section, we demonstrate the ability of the ICP
to guarantee correct registration.

C. Local Shape Matching The proposed registration algorithm is de nitely not useful
if only a subset of the data point shape P corresponds to

the model shape X or a subset of the model shape X, that is, the data point set includes a signi cant number of data
points that do not correspond to any points on the model

shape. Unfortunately, this is a trait of the majority of the shape matching algorithms that have ever been implemented. Moreover, this is a common problem in computer vision since data segmentation algorithms often misgroup one set of data points with another distinct group of points that should not be
algorithm to perform local point set matching without correspondence. Table I lists a point set with eight points that is to be matched against a set of 11 points. Fig. 4 shows the two point sets prior to registration. Fig. 5 shows the two point sets after registration by the ICP algorithm after six iterations, which took less than 1 s. The computed registration is
Translation: -48.078 6.65685 119.479 Rotation Axis: (0.0321865 0.998188 -0.0508331) Rotation Angle: 55.7188 degs RMS Error: 0.437608 mm
grouped together.
The algorithm does not pay attention to the extra points or to the ordering of the points because it always pairs a given
249
"4 fl \ 1 _,_ 1 ___' 1 '- 1-"\ 1 1 I I I, \ 1' \ f P, \ 1 T) "~\ I tr /1 I ,1
translations to achieve local matching. Compared with basic

, I r I 1 I 1' ;
point set matching, which requires the same number of points
listed in direct correspondence, we are essentially trading off

additional CPU time for local matching capability and point
/1 H
1 1
4 6
1 1
I I
/
1
J
1 I J T I I I/I I"!
, ,r I
order insensitivity.
1) Computational Issues: If one were to use a brute-force
1
1 I 1 I
tree search (testing every possible correspondence and choosing the best one) in order to
]
1 1 1 j 1 I r
/ , / /I I , j I I
nd the best match, this type
| ig
~ i
,1
I_____L___
1
I PT
of registration would require
operations of the basic
1 I
1
/ 1 1
'
1
j
1
1
1 I "1 I
f I
1/
(1
fl
least squares quaternion match. For the simple example of local matching above with Np : 8 and N,,, : 11, a brute force test would require 6 652 800 quaternion registration operations. The above registration required only six operations since the initial state was already a member of the equivalence class of the global minimum. Even with 24 initial translation states and 60 initial rotation states allowing ten iterations for
Fig. 4.
l |
Set 1 and set 2 prior to registration.

\
each combined initial state, we would require only 14 400

1 \_ 1-"'\ \
\ .._
1 1 1I
I l l
-If
Ir
I
, 1' I1
,1 , 1
\
\ I , I /
,.
-" F _, I
-_
' 11
If I !' f
i
I
/I
f , 1
1
If /I 1 I
I r " 1 "
"1 I I "
operations to provide an exhaustive and very capable matching ability. Moreover, it only takes a few minutes for these types of small point sets. The computation reduction ratio of the ICP algorithm compared with brute-force testing for this simple case is 462:1. Of course, other alternatives are possible, but we see that considering 1440 initial states is not unreasonable when the ICP algorithm is used to move from initial state to local minimum.
1
1 1- I
r
I ,
,I
I I
If
I
1 1.
1 I I I
1
1 1 I
In the African mask example in the surface section below,

we accurately registered a point set with N,,=2500 points to a point set with N,,=4200 points using 60 initial rotation
.1
1
J ___._
'
I
1'
1 - I 1 I
-. ~. \ \
"
states in less than 8 min. The amount of brute-force enumeration testing required for this case is ridiculously large; more than 17002(> 107500) operations are required. Even at 1 TeraOp/s (1012) for the age of the universe 1018 s, this exact
l
1 I 1/ " I
J1 I I
brute-force enumeration of all possible combinations would

require well over 10250 universe lifetimes!
Fig. 5.
Set 1 and set 2 after ICP registration.

TABLE I TWO 3-D POINT SETS: SET 1 IS A SUBSET OF SET 2
B. Curve Matching
22
.
u_
1' 1
..
yl
21
lT _
I2
- it l
y2 7.12 24.80 18.28
In this section, the ability of the ICP algorithm to do local free-form curve matching is demonstrated. A 3-D parametric
space curve spline was de ned as a linear combination of
43.89
42.02 42.01
-5.88
20.52
106.99 112.52
72.78 70.19 76.21
146.10
148.67
cubic B-splines and control points. A copy of it was translated

and rotated to be relatively dif cult to match. The rotated and translated curve was converted to a polyline description with 64 points. Each a:,y, z component of each point of the
25.39
4.69 17.96
113.25
112.60
147.20
148.09
44.95
44.12 48.26 Y 46.28
72.71
70.67 64.38 72.47 69.82
17.69
4.62 -5 .40 -1.16 19.81 25.00
115.15
113.59 114.58 117.65
145.95
143.59
-1.37 i 7.03 18.52
polyline was then corrupted by zero-mean Gaussian noise with a standard deviation of o : 0.1 (compared with a curve size
of 2.3 >< 2 >< 1 units). The i3o range of 0.6 units is clearly
143.85
148.32 150.00 140.00
145.00
47.06
visible compared with the size of the object. We then cut off
over half of the noisy polyline leaving a partial noisy curve shape. Fig. 6 shows the two space curves prior to registration. Fig. 7 shows the two space curves after ICP registration using
77.00
l 1
80.00
-10.0
30.00
83.00
12 initial rotation states and six initial translation states for a

total of 72 initial registration states. If the registration to move
the curve away from the original curve is post-multiplied by point on the rst set to the closest point on the other set. Only
one initial rotation state and one initial translation state were used in this example of local matching. In general, one must use an initial set of rotations combined with an initial set of the registration recovered using the ICP algorithm, we should
obtain a matrix close to the identity matrix except for the

effects of noise. The match of the partial noisy curve to the
original spline curve yielded the following registration matrix,
25O
"
"'7
"
\*
i
| y
__
'/'1'. -'. \ I '1 I / I Q I -
\ \ ml? . .\\.", ' ., I \\\ '

Q. Q \
1
I
i'\
F \\ \ .| .
i
.
\ \ \
_ _
"u '1.-/ I-~ $1.. 1-5: 'r_., ,~. . ,_ ~.:_' c_,-.,4_,-._,1 n.,'>,__-';_:,__ r,;_:,~_.__- C 4; 'r "'-2' -.~ 'F ;u -r. --r 3; us1_ -1 -.'/7f .-5/ *-,_ ;- I \ Q I1 4 iv or I) I .3, \ 1-\ f _i.r:Ii I "K. IV /' "1 I-n__/ .n , 13 -A|_{_ _.,'\/'1, _ =_'-, /,/ t.I \ I\'.. Nu m_ 9"i.' \|l| \\ \ ~ \" z \ I 11L_ ,-T! .-,-
,I \
'4;,.1" ;.~<"-#*Z=-"
\\/ v
1
"0" iii .4
'?. -~ .: .~s~'-
~\ \\
..P' _.,' _ ,, -, /Y'-1,/', ",rf it -~'A1, ,_-4. 4,--5/' :/-\ 113;" :\_-0; .'-. '/p"l/"1 >1T '_L-,.4 "51~ ~."-_'=5/r<i b;,: '~-5=."-.=!~"~35"
4,r_5' W): .r__\;1 /v -Q
.. -t-r -\ v 1'0 '1 / I;-1' , 1", __:___. ._, T_:%, ,, 4'4 t _~ 4___: | -1-_,<;;/" |\_J-\ ,. 3 1-. ~ ;1:' -. -_F'1 /Q -..r-_/,-2_- -'t~ _
-H
--c \ -.._ ',__.' _ _| _ |-
__.-
--=
__
1
\ \
I I v' \ \ .
W._/
,- ,
"Q
I
4
...
I I I,
\ \ _,-P J
r.I'_-)1
l \n\ \. Il.!.1.!sI!lif-!4. l.I!A!A!l I I ! l;\:\Il,L!.!l!L\ \ rIl 1?:7.i!il ,nIill :n'.ger :zIifl l l g./",\IM %%/
-_, u-..
Fig. 6.
Ideal and a partial noisy space curve before registration.
Fig. 8. Noisy point set and surface patch prior to registration: 250 points matched to 450 triangles.
-1
-v--_
'\
,-
\
"0- -.
T
\
0S . .i
. "315
"
\
\ \
H
Fig. 9.
I .r
\ 1
A *".="*f.~I.-{iii l:;_'f:
v
,
I, l/
,
F
a %fi" i;i %;/% jaati
Noisy point set and surface patch after registration.
via isoparametric lines indicating 450 triangles and ts in a 3 >< 3 >< 1 units box. Following the space curve example, 3-
D vector noise with a standard deviation of 0.1 units in each direction was added to the point set data to create a noisy
Fig. 7. Ideal and a partial noisy space curve after registration.
point set. The surface patch and the noisy point set are shown
which is very close to the identity matrix:

Translation: Rotation: Matrix: 0.023540 0.006925 -0.015471 0.999925 0.012227 -0.001007 0.012262 0.998647 0.050550 0.000387 0.050557 0.998721
prior to registration in Fig. 8. After running the ICP algorithm with same 24 initial rotation states for a total of about 3 min, we obtained the result shown in Fig. 9. The initial positioning
transformation multiplied by the recovered transformation for
the noisy point set yields the following approximate identity

transformation:
Translation: -0.057329 Rotation: 0.999357 Matrix: -0.033883
-0.011731
C. Surface Matching In this section, the ability of the ICP algorithm to register
0.013923 0.018430 0.034041 0.011264 0.999328 -0.013935

0.013545 0.999840
free-form surface shapes is demonstrated.

1) A Bezier Surface Patch: A simple parametric Bezier surThis demonstrates that global matching under noisy condi-
face patch was constructed for quick testing of the free-form surface matching capability of the ICP algorithm. A set of
250 randomly positioned points was evaluated in the interior of the domain of the surface patch and translated and rotated in a random manner. The points of this point set are connected
tions works quite well.

A subset of 138 noisy points was used to test the local
matching ability. The surface patch and the noisy point subset
are shown prior to registration in Fig. 10. After running the
ICP algorithm with 24 initial rotation states and six initial

translation states for a total of about 6 min, we obtained the result shown in Fig. 11. The initial positioning trans-
by lines indicating the point list sequence; they do not indicate line geometry to be matched. The surface patch is drawn
, _| . I
\ \ l
_|
J
1 1 , ll ,1 /a .
-1 _.. _-//' 5 er, ,..|._ _ I,, A , _ -'-4- 3' r7 4'.-,-.?7_'1 I ,_ _r-, r, .1 31,-,rJJ_..,._/. , 1, '1 .
.
r ._,
\ ivy, ,
\
V ,t. r , 11 a
-. _
v\ ,~ - \<\.; z\,\\\
\ P \\ i
.-
,\ \ I.
' 1, I
' |'
1 Q 0QO -III]("l;;;I;"0.:0':-:0:-ozy:.::0;J_
'01 "'0 "10 '0 0; Q0 '0 4'1; 1; 9/ 4 \\.. .:_. 1; ,;;;;$:1Z2if1:$;;.g:ggt;tg2;;:;;:j;,!fg!::;::;::zgziziaiii5E3.

lI;9 I ooooaoi ,l O0 I $0 \ \<
00,, v.,<-04.0-4 -. C 0 QQ\\\\\\\Q6; QQ D .\\\_. 41;,4- 0,0;0;0,; \ \, \ 1o;:'O::'Io:"00:-v:;33;::o'v'Q*:\\\\\::;?:~~?:\Q \ "Ia;-"I;:'lo:'l0:-';'O Q Q.0.0.0:?;//!l"Q \\\-\\_~. /l0;,'Il;,'I0;,'0-.,,P..O.OQQO.'. I If, ..\"-\ \\\\.-'\\ \, "ll""l" 0_":1;'.4. . . Q Q 1!], I Q \\\\ \\_~.: . \ Illn, ll1,, "40),;O;-QQ O /II \\.. .. o.~Q 0; I \ .. I, II 44': "I4 '0 O Q .0 0 I "0' Q Q -'0, w \\\\
,3-Z$;g;2;;2:;::32;;Z;;:gZ:.:\\\\.\-,~::iS E j 3
.:;:-. _ __-_ __ _ 4;57:;;I.;";-:-._.;2, , _, -_ - -;;\\\\\ I-o,,,'~.. .':.'~.' - ;.;-__\ \\\\\\\ '4::;a4::,::-:0:1:::;E:;:';.\\\\\:\::'_:-;\\\
II 'I \ ~ -. _ zz;qI;:O..,;$::\.Q:::5;i3i J
5.' -.4 ~;- n_ _-F hf" __-~ j-F: .:-.:,-,_n
. \\ 9/! ~-. .. -'71 ,1,-_ "X "2-_ __-1. _,,1 .:.,_,-)~, r/(I 1 1._ -',-- -7', F_._""r.'_/-.~1 ,0! _' .1-_\~--u;\;,\I'_ ,;:i9_,=_,_:_.. _7' ' 1-I ''"- d:-',- Ad'4,1
I-\_ J
, ' ~ ~*~t\ ~_ Q \ I~ - \ ,.~> ~. . . ,'I\ I > I, 1 -R It IvLn: \ \\_\:(\ I k

i}-" _.A.
C \
mi /3.". /- -' . 9&3; 3:. . , '-5. 3? * i 1
"~
1.
: WIIfff"lI''("!;"'0'o;$77:11ZZZ25/0/0:}Ir$'0f:o,.\g\:\\\\\\ Wm n 0 0' ow" ~,,, /Io I1 NW0 \ \\ 0;," II/;:'I)' ' II): fIIIzII.I I glgiff3Z;4;:41;z0:z;:I.:I:.O. \\\\ \\
I ' u!"5$"'o"o' \\\\\
J-,4 5' 41"
'$'l 41111111/l1;,''0,'I:{o,*.%'';.$551$$5;'/a,,m"o"''\\*\\\\\
Ir/I0/3'55?"1'o9"?!'5$*$'3'3'3'3'$:$$$$i/0155i$"$o"5$'''$'3:'i\\iii\ 4 I I] 0 0 oooo0 o 1 0 0 O \
III!!!
\i
pr;- '-"4 _.p'ai-!"--/ rki.4|pI '5' 7 1|-" 'bI",4','1 'f';I"4.
. _ ., _ _
*"- - -Hr
xi. ...
I'1:
.
I ___. _
Illlliiifi iiiiiiiiiiiiiiiziizgtstit;7:;;:;:;=:;:555:?:;$:$::;$;Zvt;i.ti\;t tr, 0 ;,$;;r'i/ ) "'., L'- "_:"-?'.4 3Y /lIII1';-'//I/II1,'I{l;0;'o'0'o'oti tg tgit t \\\\\-
'/film rm"-at71/1it"Iat':~:~:~:-:':.~~w.~\t\ l\\

l i$7"'i"r/'I:"It"Io"!"t"0Wtftitiity
'- - -. "'.::,'.3|;" \l\\

-'
0952-
'1-\-\_ \ ' . 11,I _ ,.
-0 _q -. I .1-:. .- \ -~~\' I -bx.-.3, -._;_'-C. \ ?, ,-., 1+ , ffitfr 1? .1 ' .._ \\\.
'*
"3Q";\>-L
*4
''
"
Y
\ '\
4
-
//J
l~2.2==" . t,t}t!"3,
W3: .
III!!_ I l!.!.!l1!L\I l.'1!4! '.!I A,!. I . 'i.It1
I'11: .-
lAm-".r'i.r/
I 151%}: I
__-
\\\\II!
Fig. 12. Model surface: Range image of mask: 8442 triangles.
Fig. 10. Noisy point subset and surface patch prior to registration: 138 points and 450 triangles.
I. -I U
_| a
it
it
1"
|
r-.
\\\.
.
2/its *1 '|l-i\ \ "

rt-
. .1 |' . ..', . ,.i. t. . t. ; . 1- t

. .
t__ -..
' L
I
\,r'..\'."~\" "
0 r\
-.;,
'\- .
\
r 1,. r !
&5E.8.%i;2:\
/ 3
-4.._
'--'.Z."'-41-PF.-"'_*-I:-_-1_\.';-_-
Fig. 11.
Noisy point subset and surface patch after registration.

_ _. _ _-_ .:.'- +5.
A -.,- __*"_'_
-=' _;~__%7T-:_~7\!>
formation multiplied by the recovered transformation for the

_.-
._ <=____ "f?TE r.-_':-
- L:- _-" \1,o.\t.:': 0 1. 0 0 Q :4:0 0,-0 0
noisy point subset yields the following approximate identity

transformation:
- "J-if/_Z -;: - 4A-J;4:
____
:_- -*1C
Fig. 13. Data: Rotated and translated point set of mask: 2546 points. Translation: -0.166476 0.159480 0.128289 Rotation: 0.990806 0.113548 0.073548 Matrix: -0.113595 0.993521 0.003545 -0.073474 -0.004842 0.997285
This matrix approximation to the identity is much less precise than the global matching case, but the data is so noisy and the shape is so featureless that we found it surprising that the registration came out as well as it did.
2) The NRCC A ican Mask: In this experiment, range data
from the National Research Council of Canadas African mask example was obtained using the commercially available Hyscan laser triangulation sensor from Hymarc, Ltd. A lowresolution 64 >< 68 gridded image was computed from the original data set for use in our experiments and is shown in
version of the measured data point set containing 2546 points is used as the data point set and is shown as scan lines in a test registration view in Fig. 13. All trial positionings of the 90-mm object, including the one shown above, converged to the correct solution in about 10 min, and all cases had a 0.59-mm RMS error. This includes six iterations worth of testing on each of 24 initial state vectors and full iteration on the best state of the six initial iterations. A side view of both the digital surface model and the measured point set as registered is shown in Fig. 14. The registration is quite accurate. A Bezier surface patch model of the mask was created to test the parametric surface matching capability on the given shape. This model is shown in in Fig. 15. The point set
Fig. 12. This data will serve as the model surface description
with 8442 triangles (4221 quadrilateral polygons). A thinned
in various rotations and translations was then registered to the parametric surface model. We had expected about a 1.2-
252
1"*
6,1,
t
up
\ \\ ..' . -r .1
i '1 fit W
4-5=";%f_-re/J
%E" T" -'=,_=*'?"-=5"
~\1t~:.~?i~:. f~\ '

i-f;i"'.
.;\\.
tn-*>.-"~'@=../AT r.~ 1
1%. 1.;-.if4? e .7
J .4_-|#"'
1'
T*"'
44- 14-
tri ( \ \-','\\i.\\
-,4 /'
-_. - -
-K\
\ \
is
In
\Q / \\ 7i-T: 1, _T.l' /
/>4
\\\\\\\
CF21? Q4; ii: e-17igfil

-
\ }YI
=
0
-~ .
A %E-- 4724/ l, i /l 1/\ mini
Fig. 14.
\-
'_. '7 \ - - /A
' ',-.
. - I
\; )1
M iF~;-_u'.
Side view of range image model and registered point set.
Fig. 15. Surface patch boundaries of a set of parametric surfaces: 97 cubic Bezier patches.
mm rms distance, but the nal solution had a 3.4-mm rms

distance. After examining the results closely, we discovered that surfaces had not been created in regions where there were
;3\_"::1:
measured points. Note the extra data points for which there was no possible surface correspondence and the evidence of misregistration in Fig. 16. Overall, this match was not bad considering the circumstances. After a post-processing step
to ignore any measured points whose point-to-surface vectors
._
; -4 I
~r.'. -_
aw
0.
5.-:i5"i* 1-" "1,
were not within a few degrees of the surface normal at the

corresponding points, the misregistration disappeared, and the
I=_=; t-'2.ttutsx Z:"?l.=~ ""'$"i'i*i*P.it{i: sasa.%s3~".=;%?,= *Isiaat 7 "~."~itt\ '- s; .t.'...: 1' p..t.rnn==_i.;.._ 'I..=;.;a.-.7-~.,.' . I "
'7"2' 4 'I1 !"-LI
it
data locked in on the surface with the expected rms distance.

Although the ICP algorithm is not designed to handle data
- Yi\ o%n.L'4'l\u|"-"-%i'"\*;-'r"f&.'-"uV3~s--", '*=r*='-5.-,=t~v"~4r-t~;<~.*:1*.%=.-.=t .:. . , -1t\,.+"@> \.

" "1-" "5 <2-E 09"! \' ' 't-'- -V -' ~ I .
points that do not correspond to the model shape, one can

conclude that a minor misgrouping of nearby points will usually have a minor effect. Another important point to keep
\I|,@*at;ia'9 i'?i~,;;&>'{ia,$\ _ \~i..tP-W. .....~';.-P-3 I"7"" ,/.'%f1b-i? if-#'im''-~""*'"""iiii:?'ilP;35t-\** ti 'i:ti ti tl&i x&\f5;?:;.;*i-,;<?: i>5E,'~ i.34~Y\-- i 1' __ ,_.._.'-an-|_|4l' ' <'_;_1I'{;l__i- Q15 - ~..i' 2-130-
?-.1 .15
in mind is that the matching algorithm does not care about

the partitioning of the composite surface model into separate surface patches.
.0. " ' ;ful11II3,',!.."F5't-I-,H,!!!\\\\l_@'_,.\&"n,|_\\\ H , 'llWHt'!.lIJ:'@'------""..wv"::.:-.-wt1 ,\. u. u1I|'mmlW\l"'.'!!.'!!.'!f--'-5----~--' .W|"l!!F"" '
1-1'
_"'
-H-'7' _.
.'.
7.. .: dft\If !4lVH'9i@@t.-.~.w~;5;f,.,,-,M,,r~,,*,';,'9;_m.. ""' =' '""iiiii5r#_"tiii$'$?""""355?.I
'1 ,I """""""-i|_.ltll-lI"'7"" I-l'l*|H'n'wl\ I, .wl"*"':;u,nug_1I F,f,"yru:lIlW_UI '\ !mn!f '
_fij E1@"',4tWY?,jg_":;;5?,1%i5A2'a%5#,g?*'
As a nal test on the mask, about 30% of the points from

the measured data point set were deleted as shown in Fig.
1./'7
~ o'er/27%~1*
'
~*' i ~ '~ 7
v i
17. The registration algorithm locked in on the solution and gave a slightly improved nns distance in less time than the full data set.
3) Terrain Data: For the nal set of experimental results, some terrain data for an area near Tucson was obtained from
Fig. 16.
Side view of parametric surface model and registered point set.
contained 13 655 points and is shown in Fig. 19. The data
the University of Arizona. Fig. 18 shows a shaded image

of the rugged terrain. The dimensions of the model surface
point set was then lifted above the model surface and rotated to be approximately orthogonal to the model set as shown
in Fig. 20. The ICP algorithm performed local matching to the model surface using 24 initial rotations and one initial translation. The registration process for these larger data sets
are 6700 >< 6840 >< 1400 units. A point set was extracted
by performing 56 planar section cuts at regular intervals
along one axis and then thinning out the data using a chord
length deviation check. An interior section of this data set was extracted such that about 60% of the surface area of the original data set was covered. The resulting data set
took about 1 hr. The results are shown in Fig. 21. The initial positioning transformation multiplied by the recovered
transformation yields the following approximate identity transformation, which demonstrates that surface matching for very
253
.
. \_:, 1/ ,1: _
ii I-*
.1
" T 1T~''=_..
i"7}'{}1\1-,|]"r._
."" .\li l"
\| M [Ir
~\\\
\\\.
C
"
f "'5 J;- \\-
1'| |i
ll
\
\
iiii W
\.
\\
,4: 5
.-._
Fig. 19.
Data: Rugged terrain near Tucson, AZ: 13 655 points.
" "' ,
4'
Fig. 17.
Subset of the data point set: 1781 points.
tr-it
Fig. 20.
Data and model of rugged terrain prior to registration.
Fig. 18.
Model surface: Rugged terrain near Tucson, AZ: 45 900 triangles.
complex surface shapes works quite well:

Translation: Rotation: Matrix: -11.330336 0.999994 -0.003309 0.000230 1.404177 0.934043 0.003308 0.000231 0.999994 0.000505 0.000505 1.000000
Fig. 21.
Data and model of rugged terrain after registration.
be incorporated if a procedure for computing the closest point

on the model shape to a given point is available. If a data shape were to come in a form other than point set form, a dense set of points on the data shape can serve as the data point set. The accelerated ICP algorithm converges to a local minimum quickly in comparison with generic nonlinear optimization methods. It is fast enough that global shape matching can be achieved using a suf ciently dense sampling of unit quaternions used as initial rotation states, and local shape matching can be achieved by combining a suf ciently dense
VII. CONCLUSIONS
The iterative closest point (ICP) algorithm is able to register
a data shape with Np points to a model shape with NI primitives. The model shape can be a point set, a set of polylines, a set of parametric curves, a set of implicit curves, a set of triangles, a set of parametric surfaces, or a set of implicit surfaces. Any other type of shape representation can
254
sampling of relevant translations. The complexity of a single

ICP iteration is O(N,,N,,) in the worst case. For N, initial translations and Nq initial rotations, the overall worst case
complexity of local matching is O(N,NqN,,N,,). The advantages of the ICP shape matching algorithm are as follows: ' It handles the full six degrees of freedom.
' It is independent of shape representation. ' The surface patch or curve segment partitioning of para-
metric or implicit entities is essentially ignored by the

matching procedure. This is important for using CAD
The disadvantages of the ICP shape matching algorithm are as follows: ' It is susceptible to gross statistical outliers unless a statistically robust method is substituted at some point (either in preprocessing or registration computation) for outlier detection. ~ The fast least squares quaternion and SVD methods are not so easily adapted to weighted least squares extensions, meaning that it is dif cult to extend the fast algorithm to allow the assignment of unequal uncertainties to points, as was done in [54]. This is not a major inconvenience
for inspection applications but would yield noticably
data in its native form without elaborate user-guided preprocessing. ' It does not require preprocessing of 3-D point data,
such as smoothing, as long as the number of statistical outliers is near zero. This is usually the case with accurate
suboptimal results for navigation laser radars and long

depth of eld triangulation systems, where uncertainty
increases signi cantly with range. ' The cost of local matching can get quite large for small
allowable occlusion percentages, e.g., 10% or less. We do
noncontact sensors used for inspection purposes.

' It does not require any derivative estimation or any local feature extraction. ' Almost all aspects of the algorithm are ideally suited
not advocate our proposed method if feature extraction

techniques will successfully solve the problem.
' The generalization to matching deformable models with

high order deformations [57] is not straightforward with-
for coarse-grain or ne-grain parallel architectures. For

large problems, even remote execution procedures and distributed le systems on networks of workstations can provide worthwhile speedup without signi cant overhead.
out, e.g., enumerating a dense set of possible deformations.
' Global matching is achieved at predictable cost based on shape complexity.

' Local matching is achieved at predictable cost based
' As an extension of the outlier rejection issue, the stated algorithm does not solve the segmentation problem, of
course. If data points from two shapes are intermixed and
on shape complexity and the percentage of allowable

occlusion.
' It can handle a reasonable amount of normally distributed vector noise, with standard deviations of up to 10% of object size demonstrated above. ' The method generalizes to rt dimensions by substituting the SVD algorithm [23] for the quaternion algorithm with the added feature that re ections are allowed. ' The method can be made statistically robust by substituting iterations of the SVD-based algorithm by Haralick er
al. [28] for the quaternion algorithm to identify outliers.
The increased computational requirements are signi cant. ' It can easily be used in conjunction with other algorithms, such as the covariance matrix alignment, which preorient the data so that fewer initial rotation states are necessary. ' Shapes with three suf ciently distinct principal moments (eigenvalues) can be globally matched at a cost of only
four initial rotation states. ' It is relatively insensitive to minor data segmentation errors as indicated by the performance of the registration of points with the African mask parametric surface model.
matched against the individual shapes, the registrations will be wrong, and the mean square distance metric will be large. This is a problem with almost all of the shape matching algorithms in the literature. ' For any given xed initial set of rotations, the global shape matching capability can be defeated even without sensor noise by constructing sea urchin or planetoid shaped objects based on the set of rotations such that correct registration cannot be guaranteed. On the other hand, for a xed set of objects and no sensor noise, one can determine an initial set of registrations in a nite amount of time such that one can guarantee registration in a nite amount of time with a suf ciently small probability of error.
' In the limit of very complicated sea urchins or perfectly
spherical planets with a single 1 um bump or in the limit of very localized matching (1% of object shape or less) on any object, the ICP algorithm degenerates to brute-force 3-D template matching. Feature extraction techniques, if possible, are preferable in such circumstances.
VIII. FUTURE DIRECTIONS
' The results of the last iteration of closest point registration can be used directly as inspection results since the distance to the closest point on a surface is computed as a byproduct. ' The accelerated ICP algorithm can achieve Newton-type
Although we have tried to present compelling results, no

method with such promise is likely to be widely accepted until
more testing has been done. We believe the algorithm for point sets, polyline sets, and triangle sets are simple enough to be
widely implemented and tested.
quadratic convergence steps at less cost than a numerical

steepest descent step. No time is spent evaluating the objective function to nd worse mean square errors off the path to the local minimum goal.
The accelerated ICP algorithm is quite ef cient, but there

is plenty of room for further computational speedup. For
example, we have not discussed the use of kd trees to ensure expected cost of O(log N,,,) rather than the worst case O(N,)
255
performance in the closest point operation. Actual testing on

parallel architectures also needs to be done. Given a large but nite amount of off-line processing, it
I [19] [11] I12] [13] I14] [15] [16] [17] [18] [19] [39]
W. Boehm, G. Farin, and J. Kahmann, A survey of curve and surface

methods in CAGD, Comput. Aided Geometric Des., vol. 1, no. 1, pp.
1-60, July 1984.

R. C. Bolles and P. Horaud, 3DPO: A three-dimensional part ori-
seems reasonable to make the following statements: For a xed

set of objects, a given level of allowable occlusion, a given
maximum possible (Gaussian) noise level, and a set of initial

registration states, it is possible to estimate the probability of registration failure by carrying out exhaustive tests. For a 1% (of the smallest shape size) noise level and a 40% maximum allowable occlusion, our tests indicate that very low probabil-
entation system, Int. J. Robotic Res., vol. 5, no. 3, pp. 3-26, Fall 1986. P. Brou, Using the Gaussian image to nd the orientation of an object,
Int. J. Robotics Res., vol. 3, no. 4, pp. 89-125, 1983.
ities of failure should be achievable with minimal extra work.

The method needs to be characterized in such detail.
A single complicated object in a set of simpler objects may

require a large set of initial registrations to handle the entire
set of objects with certainty. Unfortunately, much time is then

spent testing initial registrations with simpler objects such that one is continually homing in on the same local minimum
J. Callahan and R. Weiss, A model for describing surface shape, in Proc. Conf Comput. Vision. Part. Recogn. (San Francisco, CA), June 1985, pp. 240-247. C. H. Chen and A. C. Kak, 3DPOLY: A robot vision system for recognizing 3-D objects in low-order polynomial time, Tech. Rep. 88-48, Elect. Eng. Dept., Purdue Univ., West Lafayette, IN, 1988. R. T. Chin and C. R. Dyer, Model-based recognition in robot vision, ACM Comput. Surveys, vol. 18, no. l, pp. 67-108, Mar. 1986. C. DeBoor, A Practical Guide to Splines. New York: Springer-Verlag, 1978. C. DeBoor and K. Hollig, B-splines without divided differences,
Geometric Modeling: Algorithms and New Trends (G. Farin, Ed.), SIAM,
over and over again. With some extra bookkeeping, it may

be possible to quickly recognize a familiar local minimum well that you have already fallen into a few times and abort that iteration or to use penalty function methods to penalize walking down a path that has already been explored. The appropriate shape of the penalty functions would depend on the shape of the region of attraction in 6-D, which is dif cult to quantify and analyze.
pp. 21-28, 1987. G. Farin, Curves and Surfaces in Computer Aided Geometric Design: A Practical Guide. New York: Academic, 1990. O. D. Faugeras and M. Hebert, The representation, recognition, and locating of 3-D objects, Int. J. Robotic Res. vol. 5, no. 3, pp. 27-52, Fall 1986. T. J. Fan, Describing and recognizing 3-D objects using surface properties, Ph.D. dissertation, IRIS-237, Inst. Robotics Intell. Syst., Univ. of Southern California, Los Angeles, 1988. P. J. Flynn and A. K. Jain, CAD-based computer vision: From CAD
models to relational graphs, IEEE Trans. Patt. Anal. Machine Intell.,
It may be possible to extend the basic least squares registration solution to allow deformations (independent axis scaling
[31] [32] [2 [24] [Z5] [26] I27] [28]

29 30
and bending) of the model shapes when matching to the

data shapes. Shears and separate axis scaling can be easily
vol. 13, no. 2, pp. 114-132, 1991; see also Ph.D. Thesis, Comput. Sci. Dept., Michigan State Univ., E. Lansing, MI. E. G. Gilbert and C. P. Foo, Computing the distance between smooth objects in 3D space, RSD-TR-13-88, Univ. of Michigan, Ann Arbor, 1988. E. G. Gilbert, D. W. Johnson, S. S. Keerthi, A fast procedure for computing the distance between complex objects in 3D space, IEEE J.
Robotics Automat, vol. 4, pp. 193-203, 1988.
accomodated by allowing a general af ne transformation;

allowing even quadratic bending about the center of mass
signi cantly complicates matters. Finally, these free-form shape matching methods are likely
to be useful as part of a 3-D object recognition system.
G. H. Golub and C. F. Van Loan, Matrix Computations. Baltimore, MD: Johns Hopkins Univ. Press, 1983. W. E. L. Grimson, The combinatorics of local constraints in modelbased recognition and localization from sparse data, J. ACM, vol. 33, no. 4, pp. 658-686, 1986.
W. E. L. Grimson and T. Lozano-Prez, Model-based recognition and
ACKNOWLEDGMENT The authors wish to thank A. Morgan, R. Khetan, R. Tilove,

W. Wiitanen, R. Bartels, D. Field, W. Meyer, C. Wampler, D. Baker, R. Smith, N. Sapidis, and the reviewers for their comments.
localization from sparse range or tactile data, Int. J. Robotics Res., vol. 3, no. 3, pp. 3-35, Fall 1984. K. T. Gunnarsson and F. B. Prinz, CAD model-based localization of parts in manufacturing, IEEE Comput., vol. 20, no. 8, pp. 66-74, Aug. 1987. E. Hall, J. Tio, C. McPherson, and F. Sadjadi, Measuring curved surfaces for robot vision, Comput. vol. 15, no. 12, pp. 42-54, Dec. 1982. R. M. Haralick et al., Pose estimation from corresponding point data, inMachine Vision for Inspection and Measurement (H. Freeman, Ed.). New York: Academic, 1989.
H. Hilton, Mathematical Crystallography and the Theory of Groups of
REFERENCES
ili K. S. Arun, T. S. Huang, and S. D. Blostein, Least square tting of two
3-D point sets, IEEE Trans. Part. Anal. Machine Intell. vol. PAMI-9,
31 32 33 34
no. 5, pp. 698-700, 1987. _2i R. Bajcsy and F. Solina, Three-dimensional object representation revisited, Proc. 1st Int. Con Comput. Vision (London), June 8-11, 1989, pp. 231-240. I3: P. J. Besl, Geometric modeling and computer vision, Proc. IEEE, vol.
76, no. 8, pp. 936-958, Aug. 1988.
Movements. Oxford: Clarendon, 1963; London: Dover, 1963. B. K. P. Horn, Extended Gaussian images, Proc. IEEE, vol. 72, no. 12, pp. 1656-1678, Dec. 1984 __.._., Closed-form solution of absolute orientation using unit quaternions, J. Opt. Soc. Amer. A vol. 4, no. 4, pp. 629-642, Apr. 1987 _______, Relative orientation, A.I. Memo 994, AI Lab, Mass. Inst. Technol., Cambridge, Sept. 1987. B. K. P. Horn and J. G. Harris, Rigid body motion from range image sequences, Comput. Vision Graphics Image Processing, 1989. K. Ikeuchi, Generating an interpretation tree from a CAD model for
3-D object recognition in bin-picking tasks, Int. J. Comput. Vision, vol. 1, no. 2, pp. 145-165, 1987.
41 P. J. Besl, Active optical range imaging sensors, in Advances in Machine Vision (J. Sanz, Ed.). New York: Springer-Verlag, 1989; see also Machine Vision and Applications, vol. 1, pp. 127-152, 1989. [5] i, The free-form surface matching problem, Machine Vision for Three-Dimensional Scenes (H. Freeman, Ed.). New York: Academic, 1990. [6] P. J. Besl and R. C. Jain, Three dimensional object recognition,ACM Comput. Surveys, vol. 17, no. 1, 75-145, March 1985, [7] J. Blumenthal, Polygonizing implicit surfaces, Xerox Parc Tech. Rep. EDL-88-4, 1988. [8] B. Bhanu and C. C. Ho, CAD-based 3-D object representation for robot vision, IEEE Comput., vol. 20, no. 8, pp. 19-36, Aug. 1987.
I35] [36]
A. K. Jain and R. Hoffman, Evidence-based recognition of 3-D objects, IEEE Trans. Part. Anal. Machine Intell., vol. 10, no. 6, pp. 793-802, 1988. B. Kamgar-Parsi, J. L. Jones, and A. Rosenfeld, Registration of multiple overlapping range images: Scenes without distinctive features,
Proc. IEEE Comput. Vision Patt. Recogn. Conf (San Diego, CA), June
1989.
[37] Y. Lamdan and H. J. Wolfson, Geometric hashing: A general and

ef cient model-based recognition scheme, Proc. 2nd Int. Conf Comput. Vision (Tarpon Springs, FL), Dec. 5-8, 1988, pp. 238-251.
256
[38] I39]
[40] [41] I42: -43-
S. Z. Li, Inexact matching of 3D surfaces, VSSP-TR-3-90, Univ. of Surrey, England, 1990. P. Liang, Measurement, orientation determination, and recognition of surface shapes in range images, Cent. Robotics Syst., Univ. of California, Santa Barbara, 1987. D. Luenberger, Linear and Nonlinear Programming. Reading, MA: Addison-Wesley, 1984. A. P. Morgan, Solving Polynomial Systems Using Continuation for Engineering and Scienti c Problems. Englewood Cliffs, NJ: PrenticeHall, 1987. M. E. Mortenson, Geometric Modeling. New York: Wiley, 1985. J. Mundy et al., The PACE system, in Proc. CVPR 88 Workshop; see
also DARPA IUW.
[53] I59] [69]
B. C. Vemuri, A. Mitiche, and J. K. Aggarwal, Curvature-based representation of objects from range data, Image Vision Comput., vol. 4, no. 2, pp. 107-114, May 1986. B. C. Vemuri and J. K. Aggarwal, Representation and recognition of objects from dense range maps, IEEE Trans. Circuits Syst., vol. CAS-34, no. 11, pp. 1351-1363, Nov. 1987. A. K. C. Wong, S. W. Lu, and M. Rioux, Recognition and shape synthesis of 3-D objects based on attributed hypergraphs, IEEE Trans. Patt. Anal. Machine Intell., vol. 11, no. 3, pp. 279-290, Mar. 1.989.
-44: .45-
D. W. Murray, Model-based recognition using 3D shape alone.,

Comput. Vision Graphics Image Processing, vol. 40, pp. 250-266, 1987.
-46:
i471 I481
Press, Flannery, Teukolsky, and Vetterling, Numerical Recipes in C. Cambridge, UK: Cambridge University Press, 1988. J. Rieger, On the classi cation of views of piecewise smooth objects, Image and Vision Computing, vol. 5, no. 2, pp. 91-97, 1987. Proc. IEEE Robust Methods Workshop. Univ. of Washington, Seattle. P. Sander, On reliably inferring differential structure from 3D images, Ph.D. dissertation, Dept. of Elect. Eng., McGill Univ., Montreal,
Canada, 1988.
[491 [501
[51I
P. H. Schonemann, A generalized solution to the orthogonal procrustes probIem, Psychometrika, vol. 31, no. 1, 1966 J. T. Schwartz and M. Sharir, Identi cation of partially obscured objects in two and three dimensions by matching noisy characteristic curves,
Int. J. Robotics Res., vol. 6, no. 2, pp. 29-44, 1987.
152*
[538 I541
T. W. Sederberg, Piecewise algebraic surface patches, Comput. Aided Geometric Des., vol. 2, no. 1, pp. 53-59, 1985. B. M. Smith, IGES: A key to CAD/CAM systems integration, IEEE Comput. Graphics Applications, vol. 3, no. 8, pp. 78-83, .1983. G. Stockman, Object recognition and localization via pose clustering,
Comput. Vision Graphics Image Processing, vol. 40, pp. 361-387, 1987.
Paul J. Besl (M86) received the B.S. degree summa cum laude in physics from Princeton University, Princeton, NJ, in 1978 and the M. S. and Ph. D. degrees in electrical engineering and computer science from the University of Michigan, Ann Arbor, in 1981 and 1986, respectively. Prior to 1984, he worked for Bendix Aerospace Systems in Ann Arbor and Structural Dynamics Research Corporation, Cincinnati, OH. Since 1986, he has been a research scientist at General Motors Research Laboratories, Warren, MI, where his primary research interest is computer vision, especially the practical applications of range imaging sensors. Dr. Besl is a member of the Association for Computing Machinery. He received a Rackham Distinguished Dissertation Award for his Ph. D. work on range image understanding, which was later published in book form by Springer Verlag.
R. Szeliski, Estimating motion from sparse range data without correspondence, 2nd Int. Con Comput. Vision (Tarpon Springs, FL), Dec. 5-8, 1988, pp. 207-216. [551 G. Taubin, Algebraic nonplanar curve and surface estimation in 3space with applications to position estimation, Tech. Rep. LEMS-43, Div. Eng., Brown Univ., Providence, RI, 1988. [56] ____, About shape descriptors and shape matching, Tech. Rep. LEMS-57, Div. Eng., Brown Univ., Providence, RI, 1989. I57] D. Terzopoulos, J. Platt, A. Barr, and K. Fleischer, Elastically deformable models. Comput. Graphics, vol. 21, no. 4, pp. 205-214, July 1987.
Neil D. McKay received the B. S. degree in electrical engineering from the University of Rochester, Rochester, NY, in 1980, the M.S. degree in computer and systems engineering from Rensselaer Polytechnic Institute in 1982, and the Ph. D. degree in electrical engineering from the University of Michigan in 1985. He has been with the Computer Science Department of General Motors Research Laboratories since 1985. His interests include robot path planning and optimal control.

1992 - A Method For Registration of 3-D Shapes - ICP - OK

Uploaded by

Copyright:

Available Formats

1992 - A Method For Registration of 3-D Shapes - ICP - OK

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

1992 - A Method For Registration of 3-D Shapes - ICP - OK

Uploaded by

Copyright:

Available Formats

IEEE TRANSACFIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 14, NO.

A Method for Registration of 3-D Shapes

general, uni ed approach, which generalizes to n dimensions

and provides solutions to 1) the point-set matching problem

curve or surface derivatives, and no preprocessing of 3-D data,

The main application of the proposed method as described

and will also not be addressed. In the context of inspection

I triangle sets (faceted surfaces)

| parametric surfaces: (a:(u,v),y(u,"v),z(u,'v)).

This covers most applications that would utilize a method

shape that may correspond to a model shape, and given a

geometric shape representation, estimate the optimal rotation

the data shape minimizing the distance between the shapes

tioned above. Then, the iterative closest point (ICP) algorithm

many applications is the following question: Does a segmented

as de ned in [3] and [5] as a special case of a simple,

capabilities of the ICP registration algorithm. II. LITERATURE REVIEW

tration (pose estimation, alignment, motion estimation) of 3-D

0162-8828/92$03.00 1992 IEEE

(super)quadric models [2], [27], and 3) point sets with known

done earliest by Faugeras and his group at INRIA [18], where

alternative to reaching the local minima than our proposed

test was a simple translation along one axis: a 1-D correlation

least squares SVD registration approach, which provided a

but the block is a very simple shape. Overall, this work

statistical outliers and could theoretically be substituted for

the orthonormal matrix is strictly a positive one. A recent

the maximum eigenvalue of a 4 >< 4 matrix instead of the

relatively small. This method is much quicker than the one

the extended Gaussian image (EGI) methods allowing the

Kamgar-Parsi er al. [36] also describe a method for the

work was addressed toward piecing together terrain map data.

an approximate distance metric. Global shapes (not occluded

rotations and translations. His method forms an attributed

the registration transformation can be recovered. The method

inexact approach that allows for variability in attributes as

curve shapes, but it is unclear that the effectiveness generalizes

BESL AND MCKAY: METHOD FOR REGISTRATION OF 3-D SHAPES

any arbitrary closed-connected region in the plane. For more

III. MATHEMATICAL PRELIMINARIES

curves and surfaces, see [9], [15][17], [42], [52].

a given point on the various geometric representations listed

dtri E) = r("u.)EE -1131 4051 F'(17))-

The computations to compute the distance are not closed form

to-curve and point-to-surface distances is described below.

distance metric for an individual entity is implemented. Let

for 2' : 1, Ne. The distance between a point 13' and

the parametric entity set F is

d(r?l A) = i{jn_ir1Na} d(1Ifi.-)-

d(r, F) ,{1j_1}N} d(r1, Er)

a parametric entity is creating a simplex-based approximation

C : {F'('u.)}, one can compute a polyline L(C, 6) such that the

piecewise-linear approximation never deviates from the space

segments denoted l,, and let L :

for z' : 1,. . . , N1. The

argument value of the closest point from the line segment

more than a prespeci ed distance 6. By tagging each triangle

distance between the point 13 and the triangle t is

where u E [0, 1], o G [0,1], and w E [0,1]. The required