Ly 
IS 
j- 
l- 
/ 
>. UM 02 C 025 0 
— LA 
ISPRS Commission III, Vol.34, Part 3A , Photogrammetric Computer Vision“, Graz, 2002 
  
1. The improvement in the objective function (8), after 
some iteration step, is below a chosen value. 
2. The number of iterations exceeds a chosen constant. 
5 SIMULTANEOUS REGISTRATION WITHOUT 
CORRESPONDENCES 
Given are N clouds of 3D data points (systems 
¥1,..., XN) that partially overlap, but now we have nei- 
ther correspondences between point pairs, nor confidence 
values, as in Sec. 4. But we still assume that a good initial 
position of these N 3D point sets in a global coordinate 
system is known. Furthermore we know which of the pairs 
of systems X;, X4 actually overlap. 
In Sec. 4 each data point x; € Xj that corresponds to y; € 
X, contributes to the functional F' in Eq. (8) with the term 
Q2(x;,y;) which is quadratic in the unknowns vj; (x;) = 
(€; — En) - (ej — ex) x xi. Now we also want to set 
up such a quadratic functional for each data point x; in 
an overlapping region. We present a strategy based on a 
quadratic approximation of the squared distance function, 
cf. (Pottmann Hofer, 2002). 
If two systems 9;, X, overlap we triangulate the point 
cloud in X, and refer to the triangulated point cloud as 
T (X). First we determine those points x; € X; that ac- 
tually lie in the overlapping region, i.e., their distance to 
T (X4) is below a certain threshold. For literature on these 
thresholding techniques, see e.g. (Blais Levine, 1995, Eg- 
gert et al., 1998, Zhang, 1994). 
Now, the goal is to bring the points x; closer to the geo- 
metric shape 7'(3,), i.e., to move the points x; to lower 
levels of the squared distance function to the triangulated 
surface T' (X). In (Pottmann Hofer, 2002) it is described 
how to compute for any point x; € IR? a local quadratic 
approximant Fg x, =: F: to a smooth surface. For our ap- 
plications this has to be a nonnegative quadratic function, 
Fi(x) > 0,Vx € R?. Here we do not have a smooth sur- 
face but a triangulated point cloud T (X4). We can apply 
the results of (Pottmann Hofer, 2002), if we are able to 
locally approximate the triangulated surface by a smooth 
surface. For this we can use local quadric fits to T (X4), 
see e.g. (Yang Lee, 1999). In case that the triangulation is 
too coarse, we may first apply mesh refinement techniques 
(e.g. interpolatory subdivision) to get a sufficiently dense 
triangulation (Dyn et al., 1990, Zorin, 1997). 
Now the registration of two such overlapping point clouds 
is found by iteratively minimizing the functional 
Y Fi(xi * vj (xi), (13) 
which is quadratic in the unknowns ¢;, €;, Ck, C. If we 
do not consider points x; from one overlapping region 
(3, X.) only, but points from all overlapping regions si- 
multaneously, then we have to minimize the functional 
(13) for the unknowns e1, €, . , €N, EN. In order to fix 
certain systems Xj, one simply has to set the vectors ¢;, €i 
equal to zero. 
As a simple example for a quadratic approximant F 7 of 
the squared distance function of T (X4) at the point x;, we 
would like to present a squared point-tangent plane dis- 
tance: Let y; denote the closest point of T (Y) to xi. We 
estimate the unit normal vector n; to 7 (X) in y;, e.g. by 
computing a regression plane using the points in a certain 
neighborhood of y;. In the following we refer to the plane 
with normal vector n; through y; as the tangent plane of 
T (X4) in y;. Let d; denote the oriented distance of x; to 
this plane. If x; is already close to T (X) it is better to first 
refine the triangulation and then compute the nearest point 
and proceed as mentioned above. 
We want to minimize the squared distance of the displaced 
point x; = x; + vji (x) to the tangent plane at y;. The 
distance is given by d; 4- n;vj,(x;) and thus we get the 
functional 
Qs(xi) = Fi(x}) 7 (di  vjx(xi) : Mi)” (14) 
= (d; + (€; Em Cr) | n de det(c; — Cgk,Xi, n;))?, 
which is quadratic in the unknowns cj, €;, Cx, Ck. In 
this way we can compute a quadratic term )3(x;) for all 
points x; that have been found to lie in an overlapping re- 
gion. Hence, to simultaneously register all N point clouds 
(where at least one is fixed) we have to minimize the fol- 
lowing functional 
F=Y Qax:)- (15) 
As mentioned in Sec. 4 one has to observe the fact that the 
map x; + X; + Vjo(X;) is no Euclidean rigid body motion. 
See Sec. 4.2 for the computation of the rigid motion which 
brings x; close to tips of the vectors x; + Vjo(Xi). 
The squared distance function to the tangent plane has al- 
ready been used in several variants to the ICP algorithm, 
cf. (Bergevin et al., 1996, Chen Medioni, 1992). Our ap- 
proach of a local quadratic approximant of the squared dis- 
tance function is more general and includes the squared 
distance function to the tangent plane as a special case. 
Fig. 3 gives an example for our algorithm in the case 
N = 2, namely the registration of a point cloud to a CAD 
model. In Fig. 4 the mean squared error of the data points 
to the CAD surface after each iteration is given, both for 
our algorithm and for the standard ICP algorithm. Our al- 
gorithm converges much faster than ICP, and with respect 
to the number of iterations it is comparable to Chen and 
Medioni's method. However, we solve just a linear sys- 
tem in each iteration step, whereas Chen and Medioni in 
each step run a numerical optimization algorithm. We are 
currently working on an efficient organization of the local 
quadratic approximants in a spatial data structure in order 
to speed up the algorithm. Then, industrial inspection tasks 
are expected to come close to real time performance, since 
the computations on the approximation of the squared dis- 
tance function to the CAD model can be done in a pre- 
processing step. Further ideas for acceleration and stabi- 
lization of the algorithm can be taken from (Rusinkiewicz 
Levoy, 2001). 
A - 269