IEEE
TRANSACTIONS
ON
PATTERN
ANALYSIS
AND
MACHINE
INTELLIGENCE,
VOL.
PAMI-9,
NO.
5,
SEPTEMBER
1987
V.
Cappellini
and
A.
G.
Constantinides,
Eds.
Amsterdam,
The
Netherlands:
Elsevier,
pp.
770-775.
[3]
J.
O'Rourke,
"Motion
detection
using
Hough
techniques,"
in
Proc.
Conf.
Pattern
Recognition
and
Image
Processing,
Dallas,
TX,
1981,
p.
737.
[4]
T.
M.
Silberberg,
L.
Davis,
and
D.
Harwood,
"An
iterative
Hough
procedure
for
three-dimensional
object
recognition,"
Pattern
Rec-
ognition,
vol.
17,
no.
6,
pp.
621-629,
1984.
[5]
H.
Li,
M.
A.
Lavin,
and
R.
J.
LeMaster,
"Fast
Hough
Transform,"
in
Proc.
3rd
Workshop
Computer
Vision:
Representation
and
Con-
trol,
Bellair,
MI,
1985,
pp.
75-83.
[6]
H.
Li
and
M.
A.
Lavin,
"Fast
Hough
Transform
based
on
Bintree
data
structure,"
in
Proc.
Conf.
Computer
Vision
and
Pattern
Rec-
ognition,
Miami
Beach,
FL,
1986,
pp.
640-642.
[7]
R.
Lumia,
L.
Shapiro,
and
0.
Zuniga,
"A
new
connected
compo-
nents
algorithm
for
virtual
memory
computers,"
Comput.
Graphics
Image
Processing,
vol.
22,
pp.
287-300,
1983.
[8]
A.
Bowyer
and
J.
Woodwark,
A
programmers
geometry.
London:
Butterworth,
1983.
[9]
M.
Cohen
and
G.
T.
Toussaint,
"On
the
detection
of
structures
in
noisy
pictures,"
Pattern
Recognition,
vol.
9,
pp.
95-98,
1977.
[10]
T.
M.
Van
Veen
and
F.
C.
A.
Groen,
"Discretization
errors
in
the
Hough
Transform,"
Pattern
Recognition,
vol.
14,
pp.
137-
145,
198
1.
[11]
R.
Lumia,
"A
new
three-dimensional
connected
components
algo-
rithm,"
Comput.
Vision,
Graphics,
Image
Processing,
vol.
23,
pp.
207-217,
1983.
Least-Squares
Fitting
of
Two
3-D
Point
Sets
K.
S.
ARUN,
T.
S.
HUANG,
AND
S.
D.
BLOSTEIN
Abstract-Two
point
sets
{
pi
}
and
{
p'
};
i
=
1,
2,
9
,
N
are
re-
lated
by
p'
=
Rpi
+
T
+
Ni,
where
R
is
a
rotation
matrix,
T
a
trans-
lation
vector,
and
Ni
a
noise
vector.
Given
{
pi
}
and
{
p'
},
we
present
an
algorithm
for
finding
the
least-squares
solution
of
R
and
T,
which
is
based
on
the
singular
value
decomposition
(SVD)
of
a
3
x
3
matrix.
This
new
algorithm
is
compared
to
two
earlier
algorithms
with
respect
to
computer
time
requirements.
Index
Terms-Computer
vision,
least-squares,
motion
estimation,
quaternion,
singular
value
decomposition.
I.
INTRODUCTION
In
many
computer
vision
applications,
notably
the
estimation
of
motion
parameters
of
a
rigid
object
using
3-D
point
correspon-
dences
[1]
and
the
determination
of
the
relative
attitude
of
a
rigid
object
with
respect
to
a
reference
[2],
we
encounter
the
following
mathematical
problem.
We
are
given
two
3-D
point
sets
{
pi
};
i
-
1,
2,
,N
(here,
pi
and
p'
are
considered
as
3
x
1
column
matrices)
p>
=
Rpi
+
T
+
N,
(1)
where
R
is
a
3
x
3
rotation
matrix,
T
is
a
translation
vector
(3
x
1
column
matrix),
and
Ni
a
noise
vector.
(We
assume
that
the
ro-
tation
is
around
an
axis
passing
through
the
origin).
We
want
to
find
R
and
T
to
minimize
N
E2
=
il
p1i
(Rpi
+
T)
2.
(2)
Manuscript
received
July
2,
1986;
revised
April
9,
1987.
Recommended
for
acceptance
by
S.
W.
Zucker.
This
work
was
supported
by
the
National
Science
Foundation
under
Grant
IRI-8605400.
The
authors
are
with
the
Coordinated
Science
Laboratory,
University
of
Illinois,
Urbana,
IL
61801.
IEEE
Log
Number
8715809.
An
iterative
algorithm
for
finding
the
solution
was
described
in
Huang,
Blostein,
and
Margerum
[3];
a
noniterative
algorithm
based
on
quaternions
in
Faugeras
and
Hebert
[4].
In
this
correspondence,
we
describe
a
new
noniterative
algorithm
which
involves
the
sin-
gular
value
decomposition
(SVD)
of
a
3
x
3
matrix.
The
computer
time
requirements
of
the
three
algorithms
are
compared.
After
the
submission
of
our
correspondence,
it
was
brought
to
our
attention
that
an
algorithm
similar
to
ours
had
been
developed
independently
by
Professor
B.
K.
P.
Horn,
M.I.T.,
but
not
pub-
lished.
II.
DECOUPLING
TRANSLATION
AND
ROTATION
It
was
shown
in
[3]
that:
If
the
least-squares
solution
to
(1)
is
R
and
T,
then
{
p'
)
and
{
pi
-
Rp,
-
T
}
have
the
same
centroid,
i.
e.,
(3)
p
I=
p
,,
where
Al
N
1=
-
z
P;
-
--
Z
P;'-i~
N
i=1
N
"
t
E
"
R
N
i-
=R
l
Ni=
A
I
N
P
=
-L1
Pi.
(4)
(5)
(6)
Let
A
qi
=-pi
-
pP
(7)
(8)
We
have
2
=
q
-
Rqi
112.
=
li
Therefore,
parts:
(9)
the
original
least-squares
problems
is
reduced
to
two
(i)
Find
R
to
minimize
E2
in
(9).
(ii)
Then,
the
translation
is
found
by
t=
p'
-Rp.
(10)
In
the
next
section,
we
describe
an
algorithm
for
(i)
which
in-
volves
the
SVD
of
a
3
x
3
matrix.
III.
AN
SVD
ALGORITHM
FOR
FINDING
R
A.
Algorithm
Step
1:
From
{pi
},
{
p!
}
calculate
p,
p';
and
then
{qi
},
{
q
}.
Step
2:
Calculate
the
3
x
3
matrix
N
H-
E
qiq'
i=1
where
the
superscript
t
denotes
matrix
transposition.
Step
3:
Find
the
SVD
of
H,
H1=
UAV'.
Step
4:
Calculate
X
=
vut.
(
11)
(12)
(13)
Step
5:
Calculate,
det
(x),
the
determinant
of
X.
If
det
(x)
=
+1,
thenR
=
X.
If
det
(x)
=
-1,
the
algorithm
fails.
(This
case
usually
does
not
occur.
See
Sections
IV
and
V.)
0162-8828/87/0900-0698$01.00
©
1987
IEEE
698
Authorized licensed use limited to: Princeton University. Downloaded on January 30, 2010 at 19:45 from IEEE Xplore. Restrictions apply.
IEEE
TRANSACTIONS
ON
PATTERN
ANALYSIS
AND
MACHINE
INTELLIGENCE,
VOL.
PAMI-9,
NO.
5,
SEPTEMBER
1987
B.
Derivation
Expanding
the
right-hand
side
of
(9),
N
2=
Z
(q
-Rqi)t
(q;
-
Rqi)
N
=
Z
(q'ttq'
+
q$RtRqj
-
qtRqj
-
qR'q'R
)
i=
N
=
ZLl
(qitqi
+
qqi
-
2qtRqj).
Therefore,
minimizing
E2
is
equivalent
to
maximizing
N
F
=
qi
Rqi
=
Trace
(
Rqiq!t)
-
Trace
(RH)
(14)
where
N
H
Z
q.qit
(11)
i
=I
Lemma:
For
any
positive
definite
matrix
AA',
and
any
ortho-
normal
matrix
B,
Trace
(A4At)
2
Trace
(BAAt).
Proof
of
Lemma:
Let
ai
be
the
ith
column
of
A.
Then
Trace
(BAAt)
=
Trace
(A'BA)
=Za'(Bai).
But,
by
the
Schwarz
inequality,
at
(Bai)
s
V(ata
)(atBtBa
)
=
atai.
Hence,
Trace
(BAAt)
.
Ei
aai
=
Trace
(AAt).
Let
the
SVD
of
H
be:
H
=
UAVt
Q.E.D.
(12)
where
U
and
V
are
3
x
3
orthonormal
matrices,
and
A
is
a
3
x
3
diagonal
matrix
with
nonnegative
elements.
Now
let
X
=
VU'
(which
is
orthonormal).
(13)
We
have
XH
=
VUtUAV'
=
VAV1
(15)
which
is
symmetrical
and
positive
definite.
Therefore,
from
Lemma,
for
any
3
x
3
orthonormal
matrix
B,
Trace
(XH)
2
Trace
(BXH)
(16)
Thus,
among
all
3
x
3
orthonormal
matrices,
X
maximizes
F
of
(14).
And
if
det
(X)
=
+
1,
X
is
a
rotation,
which
is
what
we
want.
However,
if
det
(X)
=
-
1,
X
is
a
reflection,
which
is
not
what
we
want.
Fortunately,
this
degenerate
case
usually
does
not
occur.
We
shall
discuss
the
situation
in
some
detail
in
the
next
two
sec-
tions.
IV.
DEGENERACY:
NOISELESS
CASE
Assume
Ni
=
0
in
(1)
for
all
i.
Then,
obviously
there
is
a
solu-
tion
R
(which
is
a
rotation,
i.e.,
det
(R)
=
+1)
for
which
{q'
}
and
{
Rqi
}
are
congruent
and
hence
E2
=
0.
From
geometrical
considerations,
it
is
easy
to
see
that
there
are
three
possibilities.
1)
{qi
}
are
not
coplanar-Then,
the
rotation
solution
is
unique.
Furthermore,
there
is
no
reflection
X
which
can
make
E2
=
0.
Therefore,
the
SVD
algorithm
will
give
the
desired
solution.
2)
[qi
)
are
coplanar
but
not
colinear-There
is
a
unique
ro-
tation
as
well
as
a
unique
reflection
which
will
make
E2
=
0.
Therefore,
the
SVD
algorithm
may
give
either.
We
shall
see
pres-
ently
that
this
situation
can
be
easily
resolved.
3)
{qi
)
are
colinear-There
are
infinitely
many
rotations
and
reflections
which
will
make
E2
=
0.
Now
we
come
back
to
the
coplanar
case.
From
examining
the
elements
of
the
3
x
3
matrix
H,
it
can
readily
be
shown
that
the
points
{
qi
}
are
coplanar,
if
and
only
if
one
of
the
three
singular
values
of
H
is
zero.
Let
the
three
singular
values
be
XI
>
X2
>
X3
=
0.
Then
H
=
X1u1v'j
+
X2U2V2
+
0
*
U3V
(17)
where
ui
and
vi
are
columns
of
U
and
V,
respectively.
Note
that
changing
the
sign
of
U3
or
V3
will
not
change
H.
Therefore,
if
X
=
VU'
minimizes
E2,
so
does
X
'
=
V
'U'
where
V
=
[V1,
V2,
-V3].
(18)
If
X
is
a
reflection,
then
X'
is
a
rotation,
and
vice
versa.
Thus,
if
the
SVD
algorithm
gives
a
solution
X
with
det
(X)
=-1,
we
form
X'
=
V
U'
which
is
the
desired
rotation.
We
mention,
in
passing,
that
the
points
{
qi
}
are
colinear,
if
and
only
if,
two
of
the
three
singular
values
of
H
are
equal.
V.
DEGENERACY:
NOISY
CASE
If
either
{
qi
}
or
{
q!'
}
are
coplanar,
then
it
can
readily
be
shown
that
the
discussion
on
the
coplanar
case
in
Section
IV
is
still
valid,
except
of
course
now
the
minimum
of
E2
is
no
longer
zero.
Hence,
if
the
SVD
algorithm
gives
a
reflection
X
=
VU',
we
can
form
the
desired
rotation
X'
=
V
'U'.
A
special
case
of
interest
is
when
N
-
3.
Then
both
{
qi
}
and
{
q!
}
are
coplanar
point
sets.
The
situation
we
cannot
handle
is
when
the
SVD
algorithm
gives
a
solution
X
with
det
(X)
=
-1,
and
none
of
the
singular
values
of
H
is
zero.
This
means
that
neither
{
qi
}
nor
{
q'
}
are
coplanar;
yet
there
is
no
rotation
which
yields
a
smaller
E2
then
the
reflection
x.
This
can
happen
only
when
the
noise
Ni
are
very
large.
In
that
case,
the
least-squares
solution
is
probably
useless
anyway.
A
bet-
ter
approach
would
be
to
use
a
RANSAC-like
technique
(using
3
points
at
a
time)
to
combat
against
outliers
[5].
VI.
SUMMARY
OF
ALGORITHM
Using
the
procedure
of
Section
III-A,
we
obtain
X
=
VU'.
1)
If
det
(X)
=
+
1,
then
X
is
a
rotation
which
is
the
desired
solution.
2)
If
det
(X)
=
-1,
then
X
is
a
reflection.
a)
one
of
the
singular
values
(
X3,
say)
of
H
is
zero.
Then,
the
desired
rotation
is
found
by
forming
xi
=
V'Ut
where
V'
is
obtained
from
V
by
changing
the
sign
of
the
3rd
col-
umn.
b)
None
of
the
singular
values
of
H
is
zero.
Then,
conven-
tional
least-squares
solution
is
probably
not
appropriate.
We
go
to
a
RANSAC-like
technique.
VII.
COMPUTER
TIME
REQUIREMENTS
Computer
simulations
have
been
carried
out
on
a
VAX
11/780
to
compare
the
three
algorithms
(SVD,
quaternion,
iterative)
with
respect
to
time
requirements.
In
each
simulation,
a
set
of
3-D
points
{
pi
}
were
generated.
They
are
randomly
distributed
in
a
cube
of
size
6
x
6
x
6
with
center
at
(0,
0,
0).
Then
{
p!
}
were
calculated
by
rotating
{
pi
}
by
an
angle
of
750
around
an
axis
through
the
origin
with
direction
cosines
(0.6,
0.7,
0.39)
followed
by
a
trans-
lation
of
(80,
60,
70),
and
finally
by
adding
to
each
coordinate
of
the
resulting
points
Gaussian
random
noise
with
mean
zero
and
standard
deviation
0.5.
Then
the
algorithms
were
used
to
estimate
R
and
T.
The
CPU
times
used
are
listed
in
Table
I.
For
the
iterative
algorithm,
the
numbers
of
iterations
are
given
in
parentheses.
The
programs
were
written
in
C.
The
IMSL
subroutine
package
was
699
Authorized licensed use limited to: Princeton University. Downloaded on January 30, 2010 at 19:45 from IEEE Xplore. Restrictions apply.
IEEE
TRANSACTIONS
ON
PATTERN
ANALYSIS
AND
MACHINE
INTELLIGENCE,
VOL.
PAMI-9,
NO.
5,
SEPTEMBER
1987
Number
of
Point
Correspondences
TABLE
I
VAX
11/780
CPU
TIME
PER
RUN
IN
ms
Method
Used
SVD
Quaternion
Iterative
cusses
an
extension
of
the
method
to
cover
both
translational
and
ro-
tational
movements.
Index
Terms-Digital
image
processing,
fast
Fourier
transform,
im-
age
registration,
image
sequence
analysis,
motion
estimation.
3
54.6
26.6
126.8
(25)
7
41.6
32.4
108.2
(12)
11
37.0
41.0
105.2
(8)
16
39.4
45.6
94.2
(5)
20
40.4
45.2
135.0
(6)
30
44.2
48.3
111.0
(6)
used
in
finding
the
SVD
(subroutine
LSVDF)
and
in
doing
the
ei-
gen
analysis
(subroutine
EIGRS)
for
the
quaternion
method.
For
the
iterative
method,
the
initial
guess
solution
was
zero
in
all
cases.
We
observe
that
the
computer
time
requirements
of
the
SVD
and
the
quatemion
algorithms
are
comparable,
while
the
time
for
the
iterative
method
is
much
longer.
However,
in
the
iterative
method,
the
solutions
were
calculated
to
7-digit
accuracy.
If
we
can
accept
10
percent
accuracy,
then
the
number
of
iterations
are
reduced
by
a
factor
of
2
to
3.
Furthermore,
the
rate
of
convergence
can
be
increased
by
overrelaxation.
REFERENCES
[1]
S.
D.
Blostein
and
T.
S.
Huang,
"Estimating
3-D
motioni
from
range
data,"
in
Proc.
1st
Conf.
Artificial
Intelligence
Applications,
Denver,
CO,
Dec.
1984,
pp.
246-250.
[2]
D.
Cyganski
and
J.
A.
Orr,
"Applications
of
tensor
theory
to
object
recognition
and
orientation
determination,"
IEEE
Trans.
Pattern
Anal.
Machine
Intell.,
vol.
PAMI-7,
pp.
663-673,
Nov.
1985.
[3]
T.
S.
Huang,
S.
D.
Blostein,
and
E.
A.
Margerum,
"Least-squares
estimation
of
motion
parameters
from
3-D
point
correspondences,"
in
Proc.
IEEE
Conf.
Computer
Vision
and
Pattern
Recognition,
Miami
Beach,
FL,
June
24-26,
1986.
[4]
0.
D.
Faugeras
and
M.
Hebert,
'A
3-D
recognition
and
positioning
algorithm
using
geometrical
matching
between
primitive
surfaces,"
in
Proc.
Int.
Joint
Conf
Artificial
Intelligence,
Karlshrue,
West
Ger-
many,
Aug.
1983,
pp.
996-1002.
[5]
M.
Fischler
and
Bolles,
"Random
sample
consensus:
A
paradigm
for
model
fitting
with
applications
to
image
analysis
and
automated
car-
tography,"
Commun.
ACM,
vol.
24,
no.
6,
June
1981.
Registration
of
Translated
and
Rotated
Images
Using
Finite
Fourier
Transforms
E.
DE
CASTRO
AND
C.
MORANDI
Abstract-A
well-known
method
for
image
registration
is
based
on
a
conventional
correlation
between
phase-only,
or
whitened,
versions
of
the
two
images
to
be
realigned.
The
method,
covering
rigid
trans-
lational
movements,
is
characterized
by
an
outstanding
robustness
against
correlated
noise
and
disturbances,
such
as
those
encountered
with
nonuniform,
time
varying
illumination.
This
correspondence
dis-
Manuscript
received
October
11,
1984.
Recommended
for
acceptance
by
S.
W.
Zucker.
E.
De
Castro
was
with
the
Dipartimento
di
Elettronica,
Informatica
e
Sistemistica,
Universita
di
Bologna,
Viale
Risorgimento
2,
40136
Bo-
logna,
Italy.
C.
Morandi
is
with
the
Dipartiniento
di
Elettronica
ed
Automatica,
Universita
di
Ancona,
via
Brecce
Bianche.
60100
Ancona,
Italy.
IEEE
Log
Number
8715507.
I.
INTRODUCTION
Let
us
consider
a
plane
image
performing
rigid
movements
of
translation
and
rotation
within
a
rectangular
domain
C
representing
the
observed
field.
Let
the
image
be
defined
by
a
density
function
vanishing
outside
a
region
A,
whose
position
varies
with
time
t,
but
is
always
fully
contained
inside
C.
The
limitations
which
arise
when
parts
of
the
image
leave
the
observed
field
will
be
mentioned
later
on.
Let
so(x,
y)
represent
the
image
at
a
reference
time,
t
=
0,
and
st
(x,
y)
be
the
present
image,
which
is
but
a
replica
of
so(x,
y)
translated
by
(xo,
yo)
and
rotated
by
00:
st(x,
y)
=
SO(u,
v),
where
x
-
xo
=
u
cos
00
-
v
sin
00
y
-
yo
=
u
sin
0O
+
v
cos
00.
(1)
(2)
In
order
to
realign
images
so
and
st,
it
is
first
necessary
to
deter-
mine
the
translation
vector
(xo,
yo)
and
the
rotation
00
from
the
information
provided
by
so(x,
y)
and
st
(x,
y).
Several
image
reg-
istration
algorithms
are
known:
for
excellent
review
papers
see
[1]
[6].
The
present
correspondence
is
meant
as a
contribution
in
the
context
of
the
phase
correlation
technique
[7],
[8],
which
was
de-
veloped
with
reference
to
purely
translatory
displacements.
Ac-
cording
to
[7],
[8],
let
SO
(
i,
q
)
and
S,
(
(,
1
)
be
the
Fourier
trans-
forms
of
so
(x,
y)
and
s,
(x,
y).
Since
in
the
case
of
pure
translations
St
(x,
y)
=
s0(x
-
xo,
y
-
yo),
it
follows
that
St
(1
7)
=
e
-j2r(xo
+
Axo)
S
O
)
(3)
Therefore,
by
inverse
transforming
the
ratio
of
the
cross-power
spectrum
of
s,
and
so
to
its
magnitude,
StSOI/IS,So
=
exp
(
-j27r(
xo
+
r1yo)),
a
Dirac
a-distribution
centered
on
(x0,
yo)
is
obtained.
In
practice
continuous
transforms
are
replaced
by
finite
ones,
and
by
inverse
transformation
a
unity
pulse
centered
on
(x0,
vo)
is
obtained,
so
that
the
translation
is
immediately
determined.
In
this
correspondence
the
principles
of
a
generalization
of
the
phase
correlation
method
for
the
registration
of
rotated
and
trans-
lated
images
[9]
are
briefly
recalled
and
the
corresponding
numer-
ical
algorithm
is
presented.
The
effectiveness
of
the
procedure
is
then
illustrated
by
means
of
simple
experiments.
Finally,
the
edge
effects
which
arise
when
the
image
completely
fills
the
observed
field
C
are
pointed
out.
This
correspondence
is
a
part
of
research
aiming
at
the
imple-
mentation
of
an
image
stabilization
system
[10]
for
the
observation
of
the
human
retina,
in
which
the
image
viewed
by
the
TV
camera
may
be
observed
on
a
monitor
free
of
the
unavoidable
spontaneous
movements,
which
do
not
allow
any
automatic
analysis
of
dynamic
effects,
such
as
the
pulsations
frequently
observed
in
blood
vessels.
The
program
is
carried
out
in
cooperation
with
the
Ophthalmolog-
ical
Clinic
of
the
University
of
Bologna
and
the
IBM
Research
Center
in
Rome.
II.
THEORY
If
s,
(x,
y)
is
a
translated
and
rotated
replica
of
s0(x,
y),
[see
(1)],
according
to
the
Fourier
Shift
Theorem
and
the
Fourier
Ro-
tation
Theorem
[11]
their
transforms
are
related
by
St(t,
17)
=
e-j2ir"xto+,o0)
So(
t
cos
00
+
-7
sin
00,
-
sin
0o
+
-7
cos
00)
(4)
0162-8828/87/0900-0700$01.00
(g)
1987
IEEE
700
Authorized licensed use limited to: Princeton University. Downloaded on January 30, 2010 at 19:45 from IEEE Xplore. Restrictions apply.