ISPRS Commission III, Vol.34, Part 3A „Photogrammetric Computer Vision“, Graz, 2002 
  
can use the velocity vectors, and thus the squared distance 
of the displaced points is given by 
Qi(Xi, y;) = (Xi + Vjo(Xi) — y; — Vro(y:))? = (5) 
(x; — yi t (€; + €; x Xj) = (Ck + €; X ya)“ 
This term Q1(x;, y;) is a quadratic function in the un- 
knowns €;, €;, €, €; ofthe instantaneous motions applied 
to the involved systems X, and X. 
There is an alternative to Eq. (5): Instead of linearizing 
the motion of €; against X9 and the motion of X; against 
Xo, one can linearize the relative motion of X j against X. 
The velocity vector v; of a point x; € X; for this relative 
motion is given by 
Vjk(Xi) — Vjo(Xi) — Vko(xi). (6) 
Here, the distance of interest is between the point x; + 
Vjp(X;) and y; (i.e., x; is interpreted to be moving with 
system X; relative to the point y; in system X4). The 
squared distance of these two points of interest is given 
by 
Q»(Xi, yi) = (x; + Vjx(Xi) = yi? = (7) 
(x; = y, + (€; + €; x Xx) — (6 +e x x) 
The term Q»(xi, y;) is again a quadratic function in the 
unknowns €;, €;, Cx, Cg. 
We see that any pair of corresponding points gives rise to 
such a quadratic term Q4 or Q» in the involved unknown 
motion parameters. That term is a first order estimate of 
the squared distance (error) after application of the mo- 
tions. Hence, to perform the error minimization, we will 
minimize the following weighted sum 
P= > wiQ2(xi, Y;). (8) 
The weight w; is the known confidence value of the pair 
(xi, ¥;). Note that since point-to-point correspondences 
are known, both Q4 and Q» can be used. Without known 
correspondences, however, it is necessary to use the veloc- 
ity vectors v; for the relative motion of X; against X, 
see Sec. S, 
The minimization of F is mathematically simple, because 
F is a quadratic function in the unknown motion param- 
eters ¢;, ¢;. Collecting all unknowns in the vector C = 
(e1, €], Ca, Ca, , CN, Ev)", we may write F in the form 
F=CT-B-C+24-C+) wi(x_y)?. (9) 
Hence, the minimizer C of F solves the following linear 
system, 
where B is a 6 x 6N matrix and A is a IN x 6N matrix. 
Note that it is very easy to fix more than one system. Fixing 
3j just requires to set both vectors c; and €; equal to zero. 
A - 268 
4.2 Computing the actual displacements from veloci- 
ties 
In the previous subsection we have estimated the displace- 
ment vector of a point (i.e., the vector pointing from the old 
to the new position) with help of the velocity vector of an 
instantaneous motion. However, displacing points in this 
way would result in an affine mapping of the correspond- 
ing system X, and not in a rigid body motion. Although 
such affine transformations are actually used in the litera- 
ture (Bourdet Clément, 1988), we prefer to compute exact 
rigid body motions in the following way. 
It is sufficient to explain this for one moving system, which 
we denote by X, and whose instantaneous displacement is 
given by the vectors c, c. In the unlikely case that there 
is no rotational part, i.e., ¢ = 0, we are done, since then 
we have a translation with the vector €, which of course 
is a rigid body motion. Otherwise we note that the veloc- 
ity field of the instantaneous motion is uniquely associated 
with a uniform helical motion. Its axis A and pitch p can 
be computed with formula (3). The idea now is to move 
points via that helical motion approximately as far as indi- 
cated by the velocity vectors (points are now moved along 
helical paths of that motion). Note that |[e|| gives the an- 
gular velocity of the rotational part. We apply a motion to 
X which is the superposition of a rotation about the axis 
A through an angle of o, — arctan(||e||) and a translation 
parallel to A by the distance of p - a. 
A rotation through an angle of o about an axis (with unit 
direction vector a — (a;,a,,a;)) through the origin is 
known to be given by x’ = R - x with orthogonal matrix 
zin 
moo 
mii 2(b4 b» + bobs) 2(biba T bob») 
2(bi b» = bobs) moo 2(bab3 si bob1) ; 
2(b, b3 zi bob») 2(bab3 = bobı) mss 
2 2 
Mog = bg — b} + b3 — b3, mas = b% — b3 — b3 + b2, 
(11) 
where by = cos(a/2), by — a,sin(a/2), ba = 
ay sin(a/2), bs = a, sin(a/2). 
The superposition of the rotation about the axis A with 
Plücker coordinates (a, a), cf. Eq. (3), through an angle 
a and the translation parallel to A by p - « is then given by 
x — R(x — p) t (p: o)a 4 p, (12) 
where R is the matrix given above and p is an arbitrary 
point on the helical axis (e.g. p = a x a). 
4.3 Iteration and termination criteria 
With the methods from 4.1 the algorithm iteratively com- 
putes instantaneous motions of the moving systems, whose 
actual displacements are then computed as in 4.2. This it- 
erative procedure is terminated if one of the two following 
conditions is satisfied.