mporal brightness 
equence shows in- 
ng mainly in posi- 
n of the integra- 
ün yields a bet- 
flow is modeled 
neighborhood U. 
1inimization pro- 
equation system 
of the individual 
> minimization of 
hood estimation 
listributed errors 
dan (1996) show 
not hold for mo- 
motions. By re- 
ith robust statis- 
stimation of mul- 
\PPROACH 
tures within con- 
clined structures 
f spatiotemporal 
lines within a 
hborhood U can 
1e direction r = 
cular to all grey 
The direction r 
dA T 
E T3 [r1,72] 
ormulated as: 
ral gradient vec- 
ant to note that 
üivalent formula- 
y extending the 
e constraint into 
The direction r can be found by minimizing 
2 / 
Pa" / h(x — x") ((Vxtg)Tr)" dx, (6) 
which is mathematically equivalent to (2). Solving the 
quadratic terms, (6) can be written as matrix equation 
rlJr — minimum, (7) 
with the three-dimensional structure tensor 
(9x 9x) (9x 9y) | (92 91) 
J= | (9x9y) (99y) (9y 9e) | - (8) 
(9x 91) (gy 9t) (Gt 9e) 
The components of J are given by 
Jpa = (9p 9a) = / h(x — x") gpgq dX’. (9) 
Again, the spatial integration can be extended into the 
time domain for local regularization without changing 
the results of the following minimization procedure 
(Jähne, 1997). 
In order to avoid the trivial solution of (7) the con- 
straint |r|| — 1 has to be imposed on r. Solving (7) 
by the method of Lagrangian multiplicators gives the 
solution that the minimum of (7) is reached, if the 
vector r is given by the eigenvector of the tensor J to 
the minimum eigenvalue. This method is known as 
orthogonal L? approximation and can be shown to be 
mathematically equivalent to total least squares esti- 
mation. 
It is important to note that the difference between the 
least squares method of Lucas and Kanade (1981) and 
the structure tensor formulation is neither imposed by 
the formulation of the minimization nor by the exten- 
sion into the temporal domain but rather by the min- 
imization procedure. While least squares estimation 
only varies the objective function with respect to the 
two components of f, the total least squares technique 
varies all three components of the spatiotemporal vec- 
tor r under the constraint ||r||o — 1. 
4.1 Computing the structure tensor 
The implementation of the tensor components can be 
carried out very efficiently by standard image process- 
ing operators. Identifying the convolution in (9) with 
a three-dimensional spatiotemporal smoothing of the 
product of partial derivatives, each component of the 
Structure tensor can be computed as 
Jpg = B (Dy - D,), (10) 
with the 3D spatio-temporal smoothing operator B 
and the differential operator D, in the direction of the 
coordinate z,. Using a binomial operator the smooth- 
ing can be performed very. A more critical point is 
the choice of an appropriate differential operator. It 
can be shown that derivative filters optimized for a 
minimum deviation from the correct direction of the 
gradient vector reduce the error by more than an or- 
der of magnitude as compared to standard differential 
operators (Section 5). 
4.2 Eigenvalue analysis 
A standard procedure in numerical eigenvalue analysis 
is the Jacobi transformation (Press et al., 1992). The 
Jacobi method is absolutely foolproof for all real sym- 
metric matrices. This is very advantageous, because 
it does not depend on the image content. In order 
to speed up the whole eigenvalue analysis, a major 
decrease in computation time can be achieved by pre- 
selecting interesting image regions by thresholding the 
trace of the matrix, trace (J) = Jux + Jyy + Ju, for 
each point before starting the Jacobi transformation. 
4.3 Computing displacements 
Different classes of 3D spatio-temporal structures can 
be identified without explicitly solving the eigenvalue 
problem. The structure tensor contains the entire 
information on the first-order structure of the grey 
value function in a local neighborhood. By analyzing 
the rank of the matrix four different cases of spatio- 
temporal structures can be distinguished (Table 1). 
The two extreme cases of rank (J) = 0 and rank (J) = 
3 represent no apparent linear motion. 
In the case of rank (J) — 1 an already oriented im- 
age structure moves with a constant velocity. This 
is the well known aperture problem in optical flow 
computation. Only one of the three eigenvectors 
has an eigenvalue larger than zero. This eigenvec- 
tor € = (er x, ely, ee) points normal to the plane 
of constant grey value in 3D space and can be used to 
compute the normal optical flow f, : 
e 
M Eum (11) 
V e + Ei 
For rank(J) = 2 an isotropic grey value structure 
moves with a constant velocity. No aperture problem 
is present in the spatio-temporal neighborhood. The 
orientation of the 3D iso-grey-value line yields the two 
components f; and fa of the optical flow. The eigen- 
vector €, = ( €s,x3 Es,y» 6s, ) to the smallest eigenvalue 
pointing into the direction of the line corresponds to 
the initially defined vector r, and f can be computed 
as: 
re (2s. ar) - (12) 
Es,t Es,t 
707