13.472J/1.128J/2.158J/16.940J
COMPUTATIONAL GEOMETRY
Lecture 21
Dr. K. H. KoProf. N. M. Patrikalakis
Copyright c 2003 Massachusetts Institute of Technology
Contents
21 Object Matching 2
21.1 Various matching methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
21.2 Methods through localization/registration . . . . . . . . . . . . . . . . . . . . 2
21.3 Classi cation of matching methods . . . . . . . . . . . . . . . . . . . . . . . . 2
21.4 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
21.4.1 Distance metric . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
21.4.2 Distance between a point and a parametric surface . . . . . . . . . . . 4
21.4.3 Distance metric function . . . . . . . . . . . . . . . . . . . . . . . . . . 4
21.5 Matching problems : CGWOS, CPWOS, CGWS or CPWS . . . . . . . . . . . 5
21.5.1 Resolving the scaling e ects . . . . . . . . . . . . . . . . . . . . . . . . 5
21.5.2 Rotation and translation . . . . . . . . . . . . . . . . . . . . . . . . . . 5
21.6 Matching problems : IGWOS, IPWOS, IGWS or IPWS . . . . . . . . . . . . . 6
21.6.1 Iterative Closest Point (ICP) algorithm [1] for IGWOS or IPWOS . . . 6
21.6.2 ICP algorithm for scaling e ects . . . . . . . . . . . . . . . . . . . . . . 6
21.7 Matching problems : NGWOS or NPWOS . . . . . . . . . . . . . . . . . . . . 7
21.7.1 Search method [2] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
21.7.2 KH method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
21.8 Matching problems : NGWS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
21.9 Matching problems : NPWS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
21.9.1 Umbilical point method [7] . . . . . . . . . . . . . . . . . . . . . . . . . 9
21.9.2 Optimization method [7] . . . . . . . . . . . . . . . . . . . . . . . . . . 12
21.10Matching problems : o set method . . . . . . . . . . . . . . . . . . . . . . . . 13
21.10.1Distance function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
21.10.2Objective function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
21.10.3Gradient vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Bibliography 16
1
Lecture 21
Object Matching
21.1 Various matching methods
Moment method
Principal component
Contour and silhouette
New representation
Localization/registration
Miscellaneous approaches
21.2 Methods through localization/registration
A basic goal of matching through localization/registration is to nd the best rigid body trans-
formation which aligns two objects as closely as possible. The correspondence search between
two objects is a key issue in nding the best transformation for matching. Correspondence
can be established by calculating distinct features of one object and locating the same ones
on the other object. Therefore, the features have to be carefully chosen such that they are
robustly extracted and invariant with respect to various transformations. Among various fea-
tures, intrinsic di erential properties are used for matching purposes. They are independent
of parametrization and methods of representation, and only depend on the geometric shape of
the object. Moreover, they are invariant under any rigid body transformations (rotation and
translation).
21.3 Classi cation of matching methods
Two types of matching can be considered: global and partial. Simply, the global matching
is regarded as the matching for whole objects, while the partial matching is considered as
2
the matching of part of objects. Matching problems can be further categorized based on the
availability of correspondence or initial transformation information between two objects and
the application of scaling. The classi cation of matching problems is summarized in Table
21.1. In the table, acronyms are used for simpli cation as follows:
C : Correspondence information is provided.
I : Initial information on correspondence is provided.
N : No correspondence information is available.
P : Partial matching.
G : Global matching.
WOS : Without scaling.
WS : With scaling.
Global matching Partial matchingCriteria
Without scaling With scaling Without scaling With scaling
Correspondence
information CGWOS CGWS CPWOS CPWS
Initial information IGWOS IGWS IPWOS IPWS
No information NGWOS NGWS NPWOS NPWS
Table 21.1: Classi cation of matching problems
When correspondence information is provided, which is one of the types CGWOS or CP-
WOS, then a matching problem is simply reduced to calculation of the rigid body trans-
formation [3, 4]. If no correspondence is known, but a good initial approximation for the
transformation is available (IGWOS or IPWOS), then popular iterative schemes such as the
Iterative Closest Point (ICP) algorithm [1] can be employed. However, when no prior clue for
correspondence or transformation is given (NGWOS or NPWOS), the matching problem be-
comes more complicated. In this case, the solution process must provide a means to establish
such correspondence information such as in [2].
Scaling is another factor that needs to be considered separately. If a matching problem
involves scaling e ects, then direct comparison of quantitative measures cannot be used any
longer. For the global matching case, a scaling factor can be estimated by the ratio of surface
areas and applied to resolve the scaling transformation. However, when it comes to partial
matching, such area information becomes useless for the scaling factor estimation. When the
correspondence information between two objects is known (CGWS or CPWS), the scaling
factor between the objects can be easily obtained by using the ratio of Euclidean distances
between two sets of corresponding points or areas, or the ratio of the principal curvatures. If
an initial scaling value as well as a good initial approximation is provided (IGWS or IPWS),
the ICP algorithm by Besl [1] or other optimization schemes such as the quasi-Newton method
3
[12] can be employed. The problem of NGWS type can be solved by the moment method using
the principal moment of inertia and ratio of areas or volumes. No attempt, however, has been
made to solve the problems of NPWS type. For a more detailed overview, see Ko [6].
21.4 Problem statement
21.4.1 Distance metric
The Euclidean distance between two points p1 and p2 is de ned as
de(p1;p2) = jp1 p2j: (21.1)
We also de ne the minimum distance between a surface r and a point p as follows:
dsp(r;p) = minfde(p;pi);8pi 2 rg: (21.2)
21.4.2 Distance between a point and a parametric surface
Let us assume that we have a point p and a parametric surface r = r(u;v), 0 u;v 1. Then
the squared distance function is de ned as follows:
D(u;v) = jp r(u;v)j2;
= (p r(u;v)) (p r(u;v)): (21.3)
Finding the minimum distance between p and r is reduced to minimizing (21.3) within the
square 0 u;v 1. Therefore, the problem needs to be broken up into several sub-problems
which consider the behavior of the distance function at the boundary and in the interior of the
bound [10]. The sub-problems are summarized as follows: Find the minimum distances (1) in
the interior domain, (2) along the boundary curves and (3) from the corner points. Among
those minimum distances, the smallest one is chosen as the minimum distance between the
point p and the surface r. A robust calculation of the minima of the distance function (21.3)
can be achieved by the Interval Projected Polyhedron (IPP) algorithm [13, 10, 14].
21.4.3 Distance metric function
A function can be de ned using the squared distance function (21.3) to formulate a matching
problem. Suppose that we have a NURBS surface rB and an object rA represented in discrete
points or surfaces. Then, a matching problem can be stated as nding the rigid body trans-
formation (a translation vector t and a rotation matrix R) so that a global distance metric
function
= X
8p2rA
dsp(rB;( Rp + t)) (21.4)
becomes minimal, where is a scaling factor.
4
21.5 Matching problems : CGWOS, CPWOS, CGWS
or CPWS
21.5.1 Resolving the scaling e ects
Matching problems of CGWS or CPWS type involve the scaling e ects. Therefore, a scaling
factor has to be estimated so that the scaling transformation is performed, before calculating
the rigid body transformation. Since correspondence information between two objects are
available, a scaling factor can be estimated by using the ratio of the principal curvatures at
the corresponding points. After scaling has been resolved, the matching problems of CGWS
or CPWS type are treated as those of CGWOS or CPWOS type.
21.5.2 Rotation and translation
Suppose that we have two 3-tuples mi and ni(i = 1;2;3), and the correspondence information
for each point. From these points, a translation vector and a rotation matrix can be calculated.
The translation vector is easily obtained by using the centroids of each 3-tuple. The centroids
cm and cn are given by
cm = 13
3X
i=1
mi; cn = 13
3X
i=1
ni; (21.5)
and the di erence between cm and cn becomes the translation vector tT = cn cm. A
rotation matrix consists of three unknown components (the Euler angles). Since the two 3-
tuples provide nine constraints, the rotation matrix may be constructed by using some of the
constraints. But the results could be di erent if the remaining constraints are used for the
rotation matrix calculation [3]. In order to use all the constraints equally, the least squares
method may be employed [3]. The basic solution by Horn [3] is described below. Suppose that
the translation has been performed. Then what is left is to nd the rotation matrix R so that
0 =
3X
i=1
jni (Rmi)j2
=
3X
i=1
jnij2 2
3X
i=1
ni (Rmi) +
3X
i=1
jRmij2 (21.6)
is minimized. Here, D = P3i=1 ni (Rmi) has to be maximized to minimize 0. The problem
can be solved in the quaternion framework. A quaternion can be considered as a vector with
four components, i.e. a vector part in 3D and a scalar part. A rotation can be equivalently
de ned as a unit quaternion q =
h
cos( 2);sin( 2)ax;sin( 2)ay;sin( 2)az
i
which represents a rota-
tion movement around (ax;ay;az) by degree. In the quaternion framework, the problem is
reduced to the eigenvalue problem of the 4 4 matrix H obtained from the correlation matrix
M:
H =
2
66
64
s11 +s22 +s33 s23 s32 s31 s13 s12 s21
s23 s32 s11 s22 s33 s12 +s21 s31 +s13
s31 s13 s12 +s21 s22 s11 s33 s23 +s32
s12 s21 s31 +s13 s23 +s32 s33 s22 s11
3
77
75; (21.7)
5
where
M =
3X
i=1
nimTi =
2
64 s11 s12 s13s
21 s22 s23
s31 s32 s33
3
75: (21.8)
The eigenvector corresponding to the maximum positive eigenvalue is a quaternion which
minimizes the equation (21.6). An orthonormal rotation matrix R can be recovered from a
quaternion q = [q0;q1;q2;q3] by
R =
2
64 q
2
0 +q
2
1 q
2
2 q
2
3 2(q1q2 q0q3) 2(q1q3 +q0q2)
2(q1q2 +q0q3) q20 +q22 q21 q23 2(q2q3 q0q1)
2(q1q3 q0q2) 2(q2q3 +q0q1) q20 +q23 q21 q22
3
75: (21.9)
The procedure described above can also be applied to the case where more than three corre-
sponding point pairs are provided.
21.6 Matching problems : IGWOS, IPWOS, IGWS or
IPWS
21.6.1 Iterative Closest Point (ICP) algorithm [1] for IGWOS or
IPWOS
Algorithm
The point set P with Np points f~pig from the data shape and the model shape X (with
Nx supporting geometric primitives: points, lines, or triangles) are given.
The iteration is initialized by setting P0 = P, ~q0 = [1;0;0;0;0;0;0]T and k = 0. The reg-
istration 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 .
a. Compute the closest points: Yk = C(Pk;X).
b. Compute the registration: (~qk;dk) = Q(P0;Yk).
c. Apply the registration: Pk+1 = ~qk(P0).
d. Terminate the iteration when the change in mean-square error falls below a preset
threshold > 0 specifying the desired precision of the registration: jdk dk+1j< .
21.6.2 ICP algorithm for scaling e ects
When initial information on transformation is provided, the ICP method can be extended to
resolve scaling e ects in the matching problem. A scaling factor is included in the objective
function (21.4). The scaling transformation is performed at step c in the ICP algorithm. In
this case we have to provide seven initial values (three for translation, three for rotation and
one for scaling).
6
21.7 Matching problems : NGWOS or NPWOS
21.7.1 Search method [2]
Chua and Jarvis [2] developed a method to align two objects through registration assuming
no prior knowledge of correspondence between two range data sets. They use a bi-quadratic
polynomial to t data points in the local area in the least squares sense and calculate the
principal curvatures and Darboux frames. Then they construct a list of sensed data points
based on the t error. Three seed points are selected to form the list such that the area
of the triangle represented by the seed points becomes maximized to reduce any mismatch
error. Three constraints (curvature, distance and direction) are imposed to sort out possible
corresponding points out of the model data set. Then a list of transformations can be obtained
from the candidate points and an optimum transformation is selected. Various searching
algorithms are described and demonstrated in [2].
21.7.2 KH method
The overall diagram of the KH method [8] is shown in Figure 21.1. The input of the process
includes two objects and three pairs of the Gaussian and the mean curvatures at three di erent
non-collinear locations. The algorithm yields the minimum value of in the equation (21.4),
and the corresponding rotation matrix R and the translation vector t. Since no scaling e ect
is involved, we assume that a scaling factor = 1.
Step 10
Step 10 is to select three non-collinear points m1, m2 and m3 on r1 away from the boundary
of r1 where each point has di erent, Gaussian K and mean curvature H values. At mi, we
have Ki and Hi, where i = 1;2;3. Next, subdivide r2 into rational B ezier surface patches
Bj (j = 1; ;n) by inserting appropriate knots [5, 11]. Then for each rational B ezier surface
patch Bj, we express Kj and Hj in the bivariate rational Bernstein polynomial basis using
rounded interval arithmetic to formulate the problem. This allows us to use the Interval
Projected Polyhedron (IPP) algorithm [10, 13] for solving nonlinear polynomial systems. For
each pair Ki and Hi, we solve the following system of equations by the IPP technique.
Kj(u;v) = Ki K;
Hj(u;v) = Hi H; (j = 1; ;n and i = 1;2;3); (21.10)
where K and H represent the uncertainty of estimated curvatures from data points. For each
pair of Ki and Hi, a list of roots Li = (uik;vik); (k = 1; ;di) is obtained.
Step 12 (selection process)
A simple pruning search based on the Euclidean distance can be applied to the selection
process. We have three lists of candidate points, Li = (uik;vik); (k = 1; ;di) from which
one 3-tuple (n1;n2;n3)
n1 = r2(u1k;v1k); n2 = r2(u2k;v2k); n3 = r2(u3k;v3k); (21.11)
7
Find a list of possible correspondences
between two objects.
Sort out 3?tuples (selection process).
Search for the optimal 3?tuple.
10
12
14
START
END
Figure 21.1: A diagram of the KH method
is selected to satisfy the following Euclidean distance constraints simultaneously
jjm1 m2j jn1 n2jj < select;
jjm2 m3j jn2 n3jj < select;
jjm3 m1j jn3 n1jj < select; (21.12)
where select is a user-de ned tolerance.
Step 14
The correspondence information between each point mi on r1 and ni on r2 is established,
from which a list of translation vectors and a rotation matrices can be obtained. We choose a
translation vector and a rotation matrix which produces a minimum value of equation (21.4)
with s = 1, see Ko et al. [8].
21.8 Matching problems : NGWS
Since global matching is considered, scaling e ects in the matching problem of this type can
be easily resolved by using integral properties such as surface areas or volumes. After scaling
transformation, the center of mass and the principal moments are used to nd the rigid body
transformation.
8
Figure 21.2: Localized surfaces
Examples are presented to demonstrate the proposed algorithms. Solids bounded by bicu-
bic integral B-spline surface patches A and B are used. Solid A is enclosed in a rectangular
box of 25mm 23:48mm 11mm. Here, the height of solid A is 25mm. Figure 21.3 shows a
sequence of operations for matching of the two surfaces using the principal moments of inertia
of input solids. In this example, for clarity, only part of the boundary surfaces of the solids
are displayed. The smaller solid has been translated, rotated, uniformly scaled and reparam-
eterized. In Figure 21.3-(A), two boundary surfaces of the input solids are shown with their
control points. Those two surfaces have similar shape but di erent numbers of control points
and parametrization. Matching the centroids of the two solids is performed by translating the
small solid by the position di erence between the centroids, which is demonstrated in Figure
21.3-(B). The orientation of the largest principal moment of inertia of the solid A is aligned
to that of the largest one of the solid B. Similarly, the remaining two orientations are aligned
based on the values of the principal moments of inertia. After matching the orientations of the
principal moment of inertia, the two solids are aligned in their orientations as shown in Figure
21.3-(C). Figure 21.3-(D) shows that the two solids match after uniform scaling obtained from
the ratio between the volumes of the two solids, 4.651, is applied to the small solid.
21.9 Matching problems : NPWS
21.9.1 Umbilical point method [7]
The correspondence search for matching problems of NGWOS or NPWOS type only deals with
qualitative aspects. Since the!-plane is not a ected by scaling, only qualitative correspondence
can be established in the process. This implies that without a scaling factor applied, a rigid
body transformation cannot be obtained for aligning two surfaces. Therefore, a scaling factor
has to be estimated before any transformation is considered.
9
(B)
(D)(C)
(A)
Figure 21.3: Matching via integral properties
Method 1
Let us assume that we have two surfaces r1 and r2, where r1 is an approximated surface of
input data points. The overall procedure is shown in Figure 21.4.
In step 100, all generic umbilical points are located on both surfaces r1 and r2 using the
IPP algorithm [9, 10]. Non-generic umbilical points are not used in this process. If a generic
umbilical point does not exist, this procedure cannot be applied.
In step 102, the correspondence search is performed. The value ! in the complex plane is
scale-independent so that qualitative correspondences can be built from this step. Suppose that
matched pairs are denoted as mk, (k = 1; ;nk), where nk is the number of matched pairs.
Then when at least one pair is found, the next step 104 is performed. If no correspondence
is established, then the algorithm stops, implying that the umbilical point method cannot be
used in this case.
Step 104 resolves the scaling transformation. To calculate a scaling factor, the normal
curvatures are evaluated at the corresponding umbilical points on both surfaces r1 and r2.
Then the ratio between them is obtained as a scaling factor. Suppose that a surface r is scaled
with a scaling factor , denoted as r . Then the normal curvature on r is scaled to be on
r . Therefore using this relation, the scaling factor can be recovered.
In step 106, after sorting out candidate points, a rigid body transformation is estimated
by using the unit quaternion method [3]. Since the number of matched pairs is more than
three and if at least three pairs survive the selection process, the problem reduces to nding
a rigid body transformation with three known corresponding pairs. Using the methods in [3]
a rotation matrix and a translation vector can be calculated. If the matched pairs fail in the
10
START
Calculate umbilics
Correspondence pair search
using w plane
Estimate a scaling factor
Calculate a rigid body
transformation
Exist?
Transformation match the normal vectors
for each pair
Match the directions of
lines of curvature
Choose one that
minimizes the distance squared
objective function
END
mk = (w 1k , w2k)
mk < 3
Translate r 1 and
and scale r 1
No
Yes
Yes
No
100
102
104
106
108
110
112
1 and r2for r
Figure 21.4: A diagram for matching using umbilics.
11
selection process, then the algorithm goes to step 108 which deals with the matching process
with less than three matched pairs.
In step 108, the orientations of the normal vectors at the corresponding umbilical points
are aligned. First, translate the scaled surface r1 by the di erence between the positions of
the matched umbilical points. Then, align the normal vectors at the umbilical points. The
alignment of the normal vectors can be achieved by using the unit quaternion method [3]. Let
us assume that we have two normal vectors n1 and n2 at the corresponding umbilical points
for r1 and r2, respectively. The problem of matching the normal vectors can be stated as:
rotate n1 around the vector Nn = n1 n2kn1 n2k by an angle formed by n1 and n2. The angle
can be calculated by = arccos(n1;n2), see [3] for details of rotation in the quaternion frame.
In step 110, matching of lines of curvature emanating from an umbilical point is performed.
Depending on the type of the umbilical point, one (lemon) or three (star and monstar) lines
of curvature pass through the umbilical point, and each direction can be determined by the
structure of the cubic terms C(x;y) as summarized in Lecture 20. The directions can be
obtained by calculating angles of the lines of curvature with respect to the local coordinate
system at the umbilical point [9, 10]. Using the angles, vectors which indicate the directions
of lines of curvature at the umbilical point can be obtained. These vectors are calculated at
the matched umbilical points on r1 and r2. Suppose that the number of the direction vectors
on r1 is nv1 and the number of the direction vectors on r2 is nv2. Choose one vector from
r2 and align all of the vectors on r1. This process produces nv1 matched cases among which
one match is chosen that minimizes equation (21.4). This alignment is achieved by rotation
around the normal vector in the tangent plane at the matched umbilical point. Therefore, the
rotation method using the unit quaternion can be used in this process [3].
Method 2
The rigid body transformation can also be obtained by using the KH method described in
Section 4.2 after the scaling transformation is resolved. The algorithm is the same as in Figure
21.4 from step 100 through step 104. After the scaling transformation is resolved, the KH
method can be used to nd the rigid body transformation between two objects.
21.9.2 Optimization method [7]
A problem with scaling e ects can be solved with an optimization technique. Since there is
no quantitative measure that can be used to estimate a scaling value, the solution scheme has
to resort to an optimum search method which can narrow down the best estimate from the
possible set of candidate solutions.
The KH method can be treated as a function of the scaling factor which yields a value of
in equation (21.4) when a scaling factor is given. Namely, steps 10, 12 and 14 in the diagram
of Figure 21.1 are grouped as a function f such that
f = ( ;R;t); (21.13)
where is the expression given in the equation (21.4), the scaling factor, R the rotation
matrix and t the translation vector. Since the rotation matrix and the translation vector are
12
obtained from the KH method, we can reduce equation (21.13) to a function of one variable,
or f = ( ). Therefore, when is given as input, then f produces the best rigid body
transformation (a translation vector and a rotation matrix) as well as the value of the object
function de ned in equation (21.4). When the selection process fails, the tolerance select is
iteratively increased so that any triplet can be obtained. When no triplet is found, then the
function f is penalized to yield a very large value.
Using the function f( ), the problem can be formulated as a one dimensional optimization
problem to nd a scaling value which yields the minimum of f. A popular one dimension
optimization scheme, the golden section search [12] can be employed to solve it.
An initial bracket [a;b;c] of the scaling factor is provided which contains an optimum value,
and satis es f(a) >f(b) and f(c) >f(b). Suppose this bracket is I0. The golden section search
starts with I0 and continues while the size of an interval containing the optimum value is
larger than a user de ned tolerance which determines the size of the interval.
(A) (B)
Figure 21.5: A matching example of NPWS type
21.10 Matching problems : o set method
21.10.1 Distance function
The minimum distance d from each measured point to an o set of the design surface is de ned
as the distance between the points Ri and Qi minus the o set distance h, see Figure 21.6
d(Ri;Qi;h) = jRi Qi hnij (21.14)
where ni is the unit normal vector of P(u;v) at Qi given by
ni = Pu(u;v) Pv(u;v)jjP
u(u;v) Pv(u;v)jj
: (21.15)
13
Qi
Qi + nih
Ri
P(u, v)
Offset at distance h
Figure 21.6: De nition of distance function.
21.10.2 Objective function
The localization procedure minimizes the Root Mean Square Distance (RMS) between
the points ri and an o set to P(u;v) at a distance h. The RMS distance is thus given by
RMS =
sPm
i=1d2(ri;qi;h)
m (21.16)
with
d(ri;qi;h) = jri qi hnij (21.17)
where qi is the point nearest to the transformed points ri on the o set surface with
distance h from the design surface P(u;v). The position of the point qi varies together
with the motion of the point ri depending on the six motion parameters. Therefore qi
can be viewed as a function depending on the variables ; ; ;tx;ty;tz;h.
It is su cient to minimize the sum of the squared distances instead of RMS distance.
Thus the objective function for the unconstrained minimization is given by
F( ; ; ;tx;ty;tz;h) =
mX
i=1
d2(ri;qi;h): (21.18)
If we denote x = ( ; ; ;tx;ty;tz;h)T then the problem can be reduced to searching for
zeros of the gradient vector eld rF(x), i.e.
rF(x) = 0: (21.19)
The Newton-Raphson method can be used to solve the nonlinear system (21.19), namely
[H(x)] x = g(x)
x(i+1) = x(i) + x: (21.20)
14
where g(x) and [H(x)] are the gradient vector and the Hessian matrix of the objective
function F(x) given by
g(x) = rF (21.21)
[H(x)] = r2F = @
2F
@xi@xj (i;j = 1;:::;7) (21.22)
where the notation x = fxig = f ; ; ;tx;ty;tz;hg is used.
21.10.3 Gradient vector
The minimization of the objective function can be speeded up by supplying the gradient vector
of F at each iteration. Considering that
[C]tx = [C]ty = [C]tz = [C]h = 0 (21.23)
t = t = t = th = 0 (21.24)
(qi)h = 0 (21.25)
(ni)h = 0 (21.26)
h = h = h = htx = hty = htz = 0 (21.27)
hh = 1 (21.28)
the gradient vector of the objective function reduces to
rF =
2
66
66
66
66
66
66
66
66
4
F
F
F
Ftx
Fty
Ftz
Fh
3
77
77
77
77
77
77
77
77
5
=
mX
i=1
2
66
66
66
66
66
66
66
66
4
d2 (ri;qi;h)
d2 (ri;qi;h)
d2 (ri;qi;h)
d2tx(ri;qi;h)
d2ty(ri;qi;h)
d2tz(ri;qi;h)
d2h(ri;qi;h)
3
77
77
77
77
77
77
77
77
5
=
mX
i=1
2
66
66
66
66
66
66
66
66
4
2d(ri;qi;h) d (ri;qi;h)
2d(ri;qi;h) d (ri;qi;h)
2d(ri;qi;h) d (ri;qi;h)
2d(ri;qi;h) dtx(ri;qi;h)
2d(ri;qi;h) dty(ri;qi;h)
2d(ri;qi;h) dtz(ri;qi;h)
2d(ri;qi;h) dh(ri;qi;h)
3
77
77
77
77
77
77
77
77
5
(21.29)
=
mX
i=1
2
66
66
66
66
66
64
2(ri qi nih) ([C] Ri (qi) (ni) h)
2(ri qi nih) ([C] Ri (qi) (ni) h)
2(ri qi nih) ([C] Ri (qi) (ni) h)
2(ri qi nih) (ttx (qi)tx (ni)txh)
2(ri qi nih) (tty (qi)ty (ni)tyh)
2(ri qi nih) (ttz (qi)tz + (ni)tzh)
2(ri qi nih) ni
3
77
77
77
77
77
75
=
mX
i=1
2
66
66
66
66
66
64
2di [C] Ri
2di [C] Ri
2di [C] Ri
2dxi
2dyi
2dzi
2di ni
3
77
77
77
77
77
75
Unfortunately, the Hessian matrix can not be expressed explicitly. Therefore we use the quasi-
Newton’s method.
15
Bibliography
[1] P. J. Besl and N. D. McKay. A method for registration of 3D shapes. IEEE Transactions
on Pattern Analysis and Machine Intelligence, 14(2):239{256, February 1992.
[2] C. S. Chua and R. Jarvis. 3D free-form surface registration and object recognition. In-
ternational Journal of Computer Vision, 17(1):77{99, 1996.
[3] B. K. P. Horn. Closed-form solution of absolute orientation using unit quaternions. Journal
of the Optical Society of America, Series A, 4(4):629{642, 1987.
[4] B. K. P. Horn, H. M. Hilden, and S. Negahdaripour. Closed-form solution of absolute
orientation using orthonormal matrices. Journal of the Optical Society of America, Series
A, 5(7):1127{1135, 1988.
[5] J. Hoschek and D. Lasser. Fundamentals of Computer Aided Geometric Design. A.
K. Peters, Wellesley, MA, 1993. Translated by L. L. Schumaker.
[6] K. H. Ko. Algorithms for Three-Dimensional Free-Form Object Matching. PhD thesis,
Massachusetts Institute of Technology, Cambridge, MA, March 2003.
[7] K. H. Ko, T. Maekawa, and N. M. Patrikalakis. Algorithms for optimal partial matching
of free-form objects with scaling e ects. Submitted for publication, December 2002, MIT,
Department of Ocean Engineering, Cambridge, MA, 2002.
[8] K. H. Ko, T. Maekawa, and N. M. Patrikalakis. An algorithm for optimal free-form object
matching. Computer Aided Design, 35(10):913{923, September 2003.
[9] T. Maekawa, F.-E. Wolter, and N. M. Patrikalakis. Umbilics and lines of curvature for
shape interrogation. Computer Aided Geometric Design, 13(2):133{161, March 1996.
[10] N. M. Patrikalakis and T. Maekawa. Shape Interrogation for Computer Aided Design and
Manufacturing. Springer-Verlag, Heidelberg, February 2002.
[11] L. A. Piegl and W. Tiller. The NURBS Book. Springer, New York, 1995.
[12] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes
in C. Cambridge University Press, 1988.
[13] E. C. Sherbrooke and N. M. Patrikalakis. Computation of the solutions of nonlinear
polynomial systems. Computer Aided Geometric Design, 10(5):379{405, October 1993.
16
[14] J. Zhou, E. C. Sherbrooke, and N. M. Patrikalakis. Computation of stationary points of
distance functions. Engineering with Computers, 9(4):231{246, Winter 1993.
17