- Email: [email protected]
Gramians, generalized inverses, and the least-squares approximation of optical flow
JOURNAL
OF VISUAL
COMMUNICATION
AND
IMAGE
REPRESENTATION
Vol. 1, No. 1, September, pp. 3-11, 1990
Gramians, Generalized Inverses, and the Least-Squares Approximation of Optical Flow1 ROGER W. BROCKETT Pierce
Hall,
Harvard
University,
Cambridge,
Massachusetts
02138
Received May 15, 1990
significant problems in real-time motion analysis. After the early work of Koenderink and Van Doorn [7] and Longuet-Higgins and Prazdny [g], Marr and Ullman [l] began a more algorithmic line of inquiry based on the study of evolving contours and Horn and Schunck [2] initiated a study based on the derivative of the intensity. Tsai and Huang [9] and Ullman [lo] have studied motion analysis based on feature point extraction. The literature is vast; this is just a sample. There has been rather little work which attempts to unify these disparate points of view. Schalkof’f’s paper [4] discusses a distributional formulation and a least-squares method but does not pursue the consequences in depth. Our goal here is to unify and extend the common themes in the work cited above by formulating a mathematically natural and conceptually straightforward approach to the problem of determining optical flow. When it comes to applications involving evolving contours, there will be close connections with the recent work of Waxman and Ullman [ 1 l] and Waxman and Wohn [6]. The work of Bruss and Horn [ 121 and Bandyopadhyay and Aloimonos [ 131 is also relevant. We address here only the problem of determining the optical flow. That is, the problem of assigning a vector field to a subset of the image plane so as to account for the observed evolution in time. The possible ambiguities in determining what three-dimensional motion generated a given optical flow is not treated here. To use the methodology advocated here, one proceeds by first defining a set of candidate vector fields which, together with the properties of the object being studied, permit the definition of certain symmetric shape and/or reflectance matrices. Because of the orthonormalization which is implicit here, it is appropriate to call these matrices gramians. These gramian matrices reveal enough about the properties of the object in question to allow one to decide whether or not it is possible to recover the rectilinear velocities, rotational velocities, etc., and, in those cases for which it is possible, to give a robust formula for these quantities. The relationship between the linear dependence of the equations for the flow and the shape of the evolving contours is discussed in Waxman and Wohn [6] in connection with polynomial contours
This paper deals with the recovery of optical flow, that is to say, with the identification of a vector field, defined on some subset of the image plane, which accounts for the infinitesimal time evolution of the image of a particular object. Our formulation is general in that it allows for the vector field to be expressed as a linear combination of a tixed set of vector fields and it allows the measurements to include (a) the velocity of feature points, (b) the velocity normal to an evolving contour, and/or (c) the velocity tangent to an intensity gradient. The method is based on least squares and an explicit formula for the generalized inverse of a classof integral operators. It involves a gramian whose invertibility is necessary and sufficient for the identification of a unique best-fitting vector field. Various important subcases have been studied earlier and reported in the computer vision literature; the emphasis here is on the systematic development of a general bO1.
0 1990 Academic Press, Inc.
1. INTRODUCTION
In computer vision it is quite standard to have large sets of data which must somehow be reduced by extracting those features which are of greatest interest in a given context. The most naive ways of carrying out this reduction usually involve techniques which are dependent on the identification of isolated features, but these methods, being local in nature, are quite sensitive to noise. Motivated by the earlier work of Marr and Ullman [ 11, Horn and Schunck [2], Hildreth [3], Schalkoff [4], Kanatani [S], and Waxman and Wohn [6], among others, we propose and defend a precise and general methodology for setting up robust algorithms for determining optical flow based on a version of least-squares theory and a suitable (function space) form of generalized inverses. The subject of optical flow has been studied quantitatively for more than 20 years. Recent interest in the subject has been motivated by the impressive advances in computer vision hardware, which has put within reach I This research was supported in part by the Office of Naval Research under Grant NOOO14-WK0504 and by the Joint Services Electronics Program. 3
1047-3203/90 $3.00 Copyright 0 1990 by Academic Press, Inc. All rights of reproduction in any form reserved.
4
ROGER W. BROCKETT
and the vector tion of a plane. curve lying on the underlying not.
fields associated with perspective projecThey show that in this context a biquartic a plane has enough variation in it to reveal flow but that curves of lower degree do
2. PRELIMINARIES
ON LEAST SQUARES
The mathematical idea which is basic in our search for robustness, one which is standard in many areas, is nothing more than a suitable version of the least-squares method of curve fitting. Suppose that we have a set of observations depending on a parameter s:
I: II Y,(S)’ Y2b)
Y(S) =
*
Y,(S).
Suppose further that these observations are generated by a linear transformation acting on an unknown, but constant, vector
x zzs-1 P TT(s)y(s)ds. ia If S is singular then the minimizing value of x is not unique. The invertibility of the symmetric, nonnegative definite matrix S is a necessary condition for this problem to have a unique answer. (See [14, p. 711.) From a geometrical point of view, least-squares problems are solved by constructing a suitable orthogonal projection. The symmetric matrix S is diagonalized by an orthogonal transformation which sets up the orthogonal projection; i.e., it orthonormalizes the columns of T with respect to the inner product structure defined by the quadratic form in error criterion. Such matrices arise in any systematic approach to orthonormalization (e.g., the Gram-Schmitt process) and are often called gramians. Note that if the weights on the various components of the error are changed then the gramian is pre and postmultiplied by a suitable digonal matrix: S -+ DSD. Increasing the level of generality slightly now, suppose we want to minimize e = I : (T(t)x - y(t), T(t)x - y(t))& subject to constraints
Xl
on x of the form Px = 6.
x2
x=
.
Xfl
y(s)
= T(s)x,
a! 5 s 5 p.
Since this equality can be thought of as giving one equation for each s it is generally an overdetermined system and we cannot expect to solve it. However, if we seek the x which minimizes the square of the difference between the left- and right-hand sides, i.e.,
To avoid unnecessary complications, we assume that the constraints are self-consistent and expressed in a minimal way-i.e., that P is p by n and of rank p. In this case we see that for S = I” TT(t)T(t)dt a we have x = s-1
I a’ TWy(t)dt
with A being given by the solution of PS-’ I”01TT(t)y(t)dt
then there is a well-known the following lemma.
formula for x which we give in
LEMMA 1. Zf the matrix
s = I,” TT(s)T(s)ds is nonsingular
then the value of x which minimizes 7 is
+ PTA,
+ PPTh = b,
i.e., by A = (PPT)-‘(b
- PS-’ I,” TT(t)y(t)dt)
.
3. VECTOR FIELDS AND GRAMIANS
In this section we formulate which involves
and answer a question
GRAMIANS,
INVERSES,
(i) A finite collection of vector fields {Fi}y=r which have a common domain of definition in euclidean space. (ii) A distinguished subset S of the domain of definition of the Fi. (iii) A “directional” vector field T defined on S. Remark. There are two prototype situations which have been considered in the literature and which are part of the motivation for our setup. In one case S is the image of a surface with reflectance R(x, y) and T is (aR/ax, aRl8y). In the other case S is the image of a contour and T is the unit normal to the contour. In either case the set of vector fields would include the infinitesimal generators of translation in the horizontal and vertical, rotation, and possibly other vector fields as dictated by the image formation law and the surface shape. Problem. Suppose that on S we are given the dot product between T and an unknown vector field FO. Find al,a2,. . . , a,, such that the linear combination a I F1 + a2F2 + ... + a,F,, approximates F0 in the sense that
q = Is((FO- 2 aiF;,T))2dS is as small as possible. Here (F, T) denotes the ordinary dot product of the vectors F and T. The subset S will either be the union of the closure of a finite collection of open sets where dS indicates the ordinary volume measure or else S will be the union of a finite set of smooth curves where dS indicates the differential of arc length. As can be seen from the previous section, the basic ingredient in the solution of this least-squares problem is the gramian G whose 0th entry is gij = I s (Fi, T) * (4,
T)ds.
Applying the results of the previous section we can assert the following. 2. Zf the gramian G is nonsingular, then the optimal value of a = (at, a2, . . . , a,,) is unique and is given in terms of (FO, T) by LEMMA
a = G-t I s
VI, 0 [1 . . .
OPTICAL
5
FLOW
fields in the plane since this is the situation which arises in computer vision. 4.
THE VECTOR
FIELDS
OF OPTICAL
FLOW
As an object moves through space, its image flows across the image plane. If (x, y) is used to denote the horizontal and vertical coordinates, respectively, of a point on the image plane, then the velocity vector associated with the image, expressed as a function of x and y, defines a vector field on part of the image plane. There are a number of formulations of the image formation process including orthographic and projective. Rigid motions may be general three-dimensional motion or may be constrained in some way. The surfaces that rigid objects present to the camera may be planar, quadratic, or of some other form. Each combination of possibilities gives rise to vector fields on the image plane having a particular form. We illustrate using some well-known examples. 4.1. Rigid 20 Motion, Motion Plane
Zmage Plane Parallel to
This is a very simple situation. We consider a planar rigid object moving in a plane parallel to the image plane. Because the plane of motion is taken to be parallel to the image plane, all points in the motion plane are at equal depth and it is immaterial whether orthographic or projective image formation is assumed. Rigid motion in the plane is characterized by two translational velocities and one rotational velocity. In this case the velocity in the image plane is given by making a choice of CX,/3, w in the equation.
[r:J=I”
rl[:l+[;].
4.2. Rigid 30 Motion, Planar Surface
Orthographic
Projection,
In this case, it is convenient to choose a coordinate system for three-dimensional euclidean space with one axis, which we take to be the Z axis, perpendicular to the image plane. Orthogonal projection is represented by POW, Y, Z) = (X9 Y) = (x,
Yh
Wo, TW.
(Fn, T>
Zf G is not invertible, unique.
AND
then the optimal
value of a is not
Note that for the purposes of the lemma the dimensionality of the euclidean space is irrelevant; however, in the remainder of this paper we use only the case of vector
As a rigid body moves in three-dimensional space, a point vi fixed in it transforms by a formula of the type
@vi + b = sl) with 0 being an orthogonal matrix. Associated with a time-dependent euclidean transformation are a pair of three vectors, the rectilinear velocity
6
ROGER
W. BROCKETT
and the angular velocity. The angular velocity is intrinsic whereas the rectilinear velocity must be defined with respect to a choice of origin. As is well known, associated with a time-dependent orthogonal matrix o(t) is an angular velocity matrix
[ I 0
03
R=
dependent with
02
-010
on the angular velocity
0
vector (0, , 02, 03)
$ o(t) = n(t)@(t). Affine transformations form
$13
023
Y’
.
-w2
01
-w3
+[Ihbz+[I z"=1TX'11 Note that the problem of determining the full three-dimensional motion from its two-dimensional flow cannot be separated from the identification of the quantities p, q, and Z. which relate to the shape and depth. The key point here is that by differentiation of this formula one sees that the rigid motion of a plane in three-dimensional space gives rise to an affine velocity field on the image plane. To establish this fact, note that in terms of CR, as introduced above, we have
are more general. They take the
with 9 being any invertible matrix. In two dimensions there :s a three-parameter family of euclidean transformations and a six-parameter family of affine transformations. The corresponding numbers in three dimensions are 6 and 12. Orthogonal projection ignores the 2 axis, and thus the distance between points on a rigid body as projected onto the image plane will generally change as the body undergoes a rigid three-dimensional transformation. In fact, we see that POqf and POqi are related by
If we evaluate this for 0 equaling the identity matrix we get a differential version of the projection equation which involves the rotational and translational velocities
[;I=[-13 :][;I+[z:J [I +[;;]+[;]zo. X
PoOq; + POb = P&,
and thus while the coordinates of the projected image after transformation are linearly dependent on those before transformation, they also depend on the unknown Z component. If (X, Y, Z) denotes a point on the surface of an opaque object, a formula for Z in terms of X and Y is tantamount to knowing the shape of the object. Thus the above expression requires a knowledge of the shape in order to be empirically useful. If the visible surface of the object is a plane, it is described by an equation Z = pX + 4 Y + ZO. In this case rigid motion in three-dimensional space yields a (possibly singular) affine motion in the image plane. Writing
we can express this transformation
as
LPY
41
Y
The right-hand side is, indeed, an affine function of x and y. It is to be noted that the set of all possible affine functions of two variables is a six-parameter set. 4.3. Rigid 30 Motion, Planar Surface
Perspective
Projection,
As a third example of planar vector fields generated by optical flow, we consider the optical flow generated by rigid 3D motion of a planar surface as above but now with a perspective image formation law. That is to say, a point (X, Y, Z) in E3 is projected onto the image plane according to P(X, Y, Z) = (X/Z,
Y/Z) = (x, y).
GRAMIANS,
Differentiation
INVERSES,
yields v, = klZ VY
- XilZ2
= ivz - YilZ2.
Assuming that the visible surface is a plane, say Z = pX + qY + ZO, it has been shown that on the image plane v, and v?, are given by quadratic functions of x and y: v, = al + a2x + u3y + @X2 + ujxy + a6y2
AND
OPTICAL
7
FLOW
where s is arc length and IZis the unit normal to the curve. There is no need for the boundary curve to be connected. The integral may include contributions from many arcs. We propose to call this the shape grumian. Feature point tracking is not included explicitly in this formalism but can be thought of in the following way. If an object has feature points p1 , p2, . . . , pk which can be tracked, surround each of them with a small circle having radius r and center pi. Evaluate the shape gramian associated with these k circles. The limit
v, = b, + b2x + b,y + b4x2 + b5.q + bby*.
The interested reader can consult Waxman and Wohn [6] for a derivation. (Note that there are but eight free parameters.) 5.
THE DIRECTIONAL
VECTOR
FIELD
Horn and Schunck [2] and Wohn et al. [15] have reported results on using the “total” derivative of the intensity f
Z(x, y) = R,i
+ Ryj
as a means for determining ,i and j, (see, e.g., [6]). This makes sense if the object in question has a significant change in its intrinsic reflectance and if the ambient lighting is sufficiently uniform so that the convected invariance assumption is appropriate. To fit this into our setup we take S to be the image of the object in question and we take T to be the vector field T = VI
(V = gradient).
In this case the formula for the relevant gramian is go = I, (Fi, V R)(Fj,
V R)dxdy.
We propose to call this the reflectance grumiun. In the work of Marr and Ullman 111, Hildreth [3], and Waxman and Wohn [6] the problem of determining the motion of an opaque object based on the movement of its boundary is discussed. It is argued that in such a situation it is the component of the velocity which is normal to the boundary which can be determined. In this case the set S is the boundary of the object in question and T is the unit normal to the boundary. This is illustrated in Section 7. The expression for the gramian is
gu = I s (F;, n)(Fj, n)ds,
is finite and depends only on the location of the feature points. We propose to call this thefeature point grumiun. Of course it may happen that some feature points are more reliable than others. If this can be quantified, then the above definition can be modified by associating a different scaling of the radii, depending on the robustness of the associated point. These three cases are to be thought of as prototypes. There will be situations in which both reflectance and shape information are significant and situations in which feature point information needs to be combined with shape and reflectance information. In these cases the appropriate gramian may contain contributions from V Z over areas, the outward pointing normal vector over arcs, and point contributions as well. One principal motivation for formulating the problem of optical flow recovery in this generality is to be able to treat such hybrid problems. 6.
THE AFFINE
REFLECTANCE
GRAMIAN
The general principles discussed above will now be applied to some specific settings. Obviously not all possibilities can be worked out here, so we have chosen to look in some detail at a few standard problems. The problem of determining the best affme approximation to optical flow based on shape data is treated in this section. Let I, be a one-parameter family of simple closed curves in the plane. Assume that these curves are represented as a function of t (= time) and their arc length s according to
x0, 4 = $0, s) YG,
s)
=
40,
s),
with $ and 4 being piecewise differentiable functions oft and s. Let R(x, y) denote the intrinsic reflectance of an object which projects onto the image plane and boundary Ir. Let Z be the perceived intensity on the image plane.
8
ROGER
W. BROCKETT
tor fieldfr(x, y) = (YX + by + e, &4x, y) = cx + dy + J Explicitly, we assume that Z is given and we want to find a, b, c, d, e, f such that
Let F be a smooth vector field on the image plane. Consider the expression for the derivative along F $ Z(x, y) = R,Ec + Ryj
r)=
For i = fi(x, y); j, = &(x, y) we have
i=R,fi
ARG=I\
is minimized. Define a 6 x 6 symmetric matrix AZ?G by the expression (the integration is over the interior of F)
R2x2 x
R:XY
R,Ryx2
R,R,xy
R:x
R,R,x
R;xy
R:y2
R,R,xy
R,R,y2
R:x
RxRyx
R,Ryx2
R,R,xy
R;XY
R?
R,R,x
R,R,x
RxR,xy
R,R,y=
R;XY
R;
RxR,y
RyRyy
R;x
R:x
RxRyx
R,R,x
R,R,x
R;x
THEOREM 1. Zf ARG is nonsingular, (a, 6, c, d, e, f) which minimizes 11 is
then the value of
zxx LY ZYX
ZYY 1.X
the minimizing
R,R,
R::
R,Ry
Rf
dxdy.
lem if ARG is nonsingular. The nonsingularity of ARG and, in fact, the size of the smallest eigenvalue of ARG, influences the robustness of this solution. We will call ARG the afine rejiectance gramian associated with Z and R. Later on we will discuss in qualitative terms the nonsingularity of ARG.
7. THE AFPINE SHAPE GRAMIAN Suppose that Z(x, y) is a function which is one on R and zero outside R. Such a case corresponds to a silhouette moving across a plain background. All the information on the motion comes from the moving edge. Let (n,. . ,., n,),, be the components of the outward pointing normal to F at s. This leads to a formula involving the symmetric matrix
value of (a, b, c, d,
This result is an immediate consequence of the results cited above. It provides a complete solution of the prob-
1
n2x2 Y
AGS = 1
R:
R&Y
-
k, y)drdy.
ZY Zf ARG is singular, e, f) is not unique.
+ by + e)
+Ryf2.
The first problem we consider is that of associating to of F by an affine vec-
11,
R (i - R&x
- R,(cx + dy + f))2 dxdy
R the least-squares approximation
= (ARG)-’
II
+y
nxnyx
2
nxnyw 2
4xy
n?y2
nxnyxy
wyx2
nay-y
n2x2
d-v
nxnyxY
nxnyY
nZxy
n2y2
nix
nix
nay
vy
n2w
+
wyx
wyx
vyY
nfx
nxnyx
n,2x
w,x
nay
+
w,y
4
vy
ds
4y 2
nY
:
wy 4
2
I
GRAMIANS,
INVERSES,
which we call the u~@e shape gramian. The leastsquares fit thus results in a formula for the motion parameters . THEOREM 2. Zf the ASG is nonsingular then the value of (a, b, c, d, e, f) which minimizes 7) is
AND
OPTICAL
9
FLOW
From the above congruent formula, we see that the shape gramian of a parallelogram is also nonsingular since affine transformation will take a square into a parallelogram. On the other hand, the ASG of an ellipse is singular since that of a circle is singular and an affine transformation can be found to convert an ellipse into a circle.
w nYX
= (ASG)-’ 1,.
w
8.
v,(s)ds.
TWO-DIMENSIONAL RIGID MOTION THE IMAGE PLANE
PARALLEL
TO
w nY n, I-
Zf ASG is singular the minimization not unique.
of (a, b, c, d, e, f) is
This formula is the appropriate one when only the edges or silhouette of an object can be observed and no shading information is available. Since the ASG can be difficult to evaluate, it is useful to discuss how it depends on coordinate transformations. The key point is that it transforms by a congruence which depends on the underlying transformation of coordinates in a simple way. THEOREM 3. Zf ASG is the ufJine shape grumiun with respect to a coordinate system (x, y) and if (x’, y’) is related to (x, y) by
[t:]=[:iI [:1+[:19 then the ufJine shape grumian with respect to (x’, y’) is given by PTSP, with P being nonsingular if and only if the mapping from (x, y) to (x’, y”) is invertible. It is possible to express P quite briefly in terms oftensoral transformation laws which make the naturality of this construction more explicit. In order to fix ideas we evaluate the affine shape gramian for a square, chasing a coordinate system centered at the center of the square and aligned with the edges. If the square is two units on a side the ASG is
This section is a digression devoted to the study of the special case of two-dimensional euclidean motion. This is of interest in its own right and as a simplification of the affine case treated above. Consider a curve in two-dimensional euclidean space which is given by {x, y[ f (x, y) = O}. Pick an orientation for the curve and let s denote arc length along the curve. Because the group of rigid motions in the plane is threedimensional and no projection is involved, we associate with such a curve a 3 x 3 symmetric matrix which is the appropriate one for extracting the translational and rotational velocities associated with rigid motion. The process involves two steps. First, we introduce a construction which involves an arbitrary choice origin. In step two we will eliminate this dependence on how the curve is displaced from the origin by associating to the curve a standard choice of origin. If the curve undergoes a rigid motion caused by an xvelocity ,, a y-velocity y, and a rotation 8 about the origin, then at a point (x, y) on the curve, the curve moves in its normal direction with a velocity
l L”]*([:]+[1,lq vw=l@p fy That is, the normal velocity is 1
EXAMPLE.
400000 040000
1
000040 00000Y
v(s) = *
+ f;
.((fc 3f, 3 Yfx
- xfyh t-k j, f%>.
This equation is of the form 4s) = (P(S), i).
1
1
If we regard the various values of v(s) as given then we can ask how best to recover the values of the parameters ,, y, 6 which characterize the rigid motion. Assuming we have scaled x and y so that it is reasonable to minimize the sum of the squares of the errors in ,, y, and 0, then the results of Section 2 apply. In this case we have
10
ROGER
W. BROCKETT
not determine i, j, 8, uniquely. It is not difficult to see that S(F) is singular if and only if I’ has constant curvature. That is, I must be a line segment or an arc of a circle. The choice of coordinate frame for this calculation has, until now, been arbitrary. We should, thus, ask how S(F) changes if the coordinate system is changed? Say we move the origin and rotate coordinates. Introduce
with
[1 PI(S)
ESG
=
Ir
p2W
[PI@),
p2bh
mWld~.
P3b)
Y’
We call ESG the two-dimensional euclidean shape gramian of the curve relative to the given frame F. If it is singular then the normal component of the velocity does
cos I30 S(F) =
If we partition
::sHBu] [:I+[I. [ I=[::iYio X’
Going back to the definition we see that this gives a new s:
-sin 80 cos 130 y0
y.
S as
rd
s = I”:,,1
a
Ls then it becomes
S=
1
=r @TOT
L(OTa + Om)
(OTa + Om)T r - 2md + aTTa
Thus the determinant and trace of Tare intrinsic. Moreover, since there is always an a such that Ta + m = 0 we can always choose an origin such that the ESG is diagonal. It is not too difficult to see that except for curves whose normal never changes this choice is unique. For straight line segments this constrains the origin to the locus of points defined by the perpendicular bisector of the segment. We refer to this choice of origin as the centroid of the curve and refer to the axes which diagonalize T as the symmetry axes of the curve. EXAMPLE 1. Consider a straight line of length r. If the origin of the coordinate system is taken to be the midpoint of the line and if the line has slope CYthenfr(s) = xl-, f2(s) = ax/V?%?, and s goes from -r/2 to r/2. In this case
1
x0 T
sin 80 ~0’
.
a
s
1
-03
-as
s2
1 ds
1 a2 ff 1 s=r[(y 1 [
a2/a2 + 1 arta + 1
0
(r/a2 + 1
l/O? + 1
0
0
0
r2/12
.
Note that in this case S is singular. If we move the origin to the endpoint of the line then S becomes
a2 + 1
1
1
r/2
r/2
r2/3
,
which is, of course, also singular. EXAMPLE 2. Consider an equilateral triangle with coordinate system at the centroid of the triangle. In terms of the perimeter p, S is given by
LO
0
p/271
GRAMIANS,
INVERSES,
EXAMPLE 3. A regular n-gon of perimeter p has a centroid located at its center and has gramian shape
s=
-p/2
0
0
0
p/2
0
-Q
0
p3/3n2
I-
We note that the ESG of a plane curve expressed with respect to the natural coordinate system is invariant with respect to rigid motion and thus can be expressed in terms of the arc length and the curvature, alone. For example, the trace of T, defined above, is the arc length.
9. CONCLUSIONS
Relative motion between the objects being viewed and the camera yield a changing pattern on the image plane. Depending on the situation, the piecewise smooth flow can be represented as a member of a parametrized family of vector fields of a specific type. Solving the inverse problem, i.e., the problem of finding the relative motion in three-dimensional space given the changing patterns in the image plane, often involves the computation of a representation of the two-dimensional flow field in the image plane. We have shown here that there is a natural common framework for this process. It involves a generalized inverse of an integral operator and through the mathematical formulation of this process we are lead to necessary and sufficient conditions for the determination of the optical flow.
AND OPTICAL
II
FLOW
REFERENCES 1. D. Marr and S. Ullman, Directional selectivity and its use in early visual processing, Proc. Roy. Sot. London Ser. B. 211, 1981, 151180. 2. B. K. Horn and B. G. Schunck, Determining optical flow, ArtiJicial Intelligence 17, 1981, 185-208. 3. E. C. Hildreth, Computations underlying the measurement of visual motion, Art$cial Intelligence 23, 1984, 309-354. 4. R. J. Schalkoff, Image Motion Analysis Using the Concept of Weak Solutions to Distributed Parameter Systems, CVPR, 1983. 5. K. Kanatani, Structure from motion without correspondence: General principle, Proc. IJCHA, 1985, 886-888. 6. A. M. Waxman and K. Wohn, Contour evaluation, neighborhood deformation and global image flow: Planar surfaces in motion, Internat. J. Robotics Res. 4, 1985, 95-108. 7. J. J. Koenderink and A. J. Van Doom, Local structure of movement parallax of the plane, J. Opt. Sot. Amer. 66, 1976, 717-723. 8. H. C. Longuet-Higgins and K. Prazdny, The interpretation of a moving retinal image, Proc. Roy. Sot. London Ser. B. 208, 1980, 385-397. 9. R. Y. Tsai and T. S. Huang, Uniqueness and estimation of threedimensional motion parameters of rigid objects with curved surfaces, in Proceedings, PRIP 1982, Las Vegas, NV, pp. 112-118. IO. S. Ullman, The Interpretation of Visual Motion, MIT Press, Cambridge, MA, 1979. Il. A. M. Waxman and S. Ullman, Surface structure and 3D motion from image flow: A kinematic analysis, Znternnt. J. Robotics Res. 4, 1985, 72-94. 12. A. R. Bruss and B. K. P. Horn, Passive navigation, Comput. Vision, Graphics Image Process. 2, 1983, 3-20. 13. A. Bandypopadhyay and J. Aloimonos, Perception of rigid motion from spatio-temporal derivatives of optical flow, Unpublished manuscript . 14. R. W. Brockett, Finite Dimensional Linear Systems, Wiley, New York, 1970. 15. K. Wohn, L. S. Davis, and P. Thrift, Motion estimation Based on multiple local constraints and nonlinear smoothing, Pattern Recognition 16, 6, 1983, 563-570. 16. E. J. DeIp and P. H. Eichel, Motion parameter estimation: Developing the scale change parameter, Unpublished manuscript.
OF VISUAL
COMMUNICATION
AND
IMAGE
REPRESENTATION
Vol. 1, No. 1, September, pp. 3-11, 1990
Gramians, Generalized Inverses, and the Least-Squares Approximation of Optical Flow1 ROGER W. BROCKETT Pierce
Hall,
Harvard
University,
Cambridge,
Massachusetts
02138
Received May 15, 1990
significant problems in real-time motion analysis. After the early work of Koenderink and Van Doorn [7] and Longuet-Higgins and Prazdny [g], Marr and Ullman [l] began a more algorithmic line of inquiry based on the study of evolving contours and Horn and Schunck [2] initiated a study based on the derivative of the intensity. Tsai and Huang [9] and Ullman [lo] have studied motion analysis based on feature point extraction. The literature is vast; this is just a sample. There has been rather little work which attempts to unify these disparate points of view. Schalkof’f’s paper [4] discusses a distributional formulation and a least-squares method but does not pursue the consequences in depth. Our goal here is to unify and extend the common themes in the work cited above by formulating a mathematically natural and conceptually straightforward approach to the problem of determining optical flow. When it comes to applications involving evolving contours, there will be close connections with the recent work of Waxman and Ullman [ 1 l] and Waxman and Wohn [6]. The work of Bruss and Horn [ 121 and Bandyopadhyay and Aloimonos [ 131 is also relevant. We address here only the problem of determining the optical flow. That is, the problem of assigning a vector field to a subset of the image plane so as to account for the observed evolution in time. The possible ambiguities in determining what three-dimensional motion generated a given optical flow is not treated here. To use the methodology advocated here, one proceeds by first defining a set of candidate vector fields which, together with the properties of the object being studied, permit the definition of certain symmetric shape and/or reflectance matrices. Because of the orthonormalization which is implicit here, it is appropriate to call these matrices gramians. These gramian matrices reveal enough about the properties of the object in question to allow one to decide whether or not it is possible to recover the rectilinear velocities, rotational velocities, etc., and, in those cases for which it is possible, to give a robust formula for these quantities. The relationship between the linear dependence of the equations for the flow and the shape of the evolving contours is discussed in Waxman and Wohn [6] in connection with polynomial contours
This paper deals with the recovery of optical flow, that is to say, with the identification of a vector field, defined on some subset of the image plane, which accounts for the infinitesimal time evolution of the image of a particular object. Our formulation is general in that it allows for the vector field to be expressed as a linear combination of a tixed set of vector fields and it allows the measurements to include (a) the velocity of feature points, (b) the velocity normal to an evolving contour, and/or (c) the velocity tangent to an intensity gradient. The method is based on least squares and an explicit formula for the generalized inverse of a classof integral operators. It involves a gramian whose invertibility is necessary and sufficient for the identification of a unique best-fitting vector field. Various important subcases have been studied earlier and reported in the computer vision literature; the emphasis here is on the systematic development of a general bO1.
0 1990 Academic Press, Inc.
1. INTRODUCTION
In computer vision it is quite standard to have large sets of data which must somehow be reduced by extracting those features which are of greatest interest in a given context. The most naive ways of carrying out this reduction usually involve techniques which are dependent on the identification of isolated features, but these methods, being local in nature, are quite sensitive to noise. Motivated by the earlier work of Marr and Ullman [ 11, Horn and Schunck [2], Hildreth [3], Schalkoff [4], Kanatani [S], and Waxman and Wohn [6], among others, we propose and defend a precise and general methodology for setting up robust algorithms for determining optical flow based on a version of least-squares theory and a suitable (function space) form of generalized inverses. The subject of optical flow has been studied quantitatively for more than 20 years. Recent interest in the subject has been motivated by the impressive advances in computer vision hardware, which has put within reach I This research was supported in part by the Office of Naval Research under Grant NOOO14-WK0504 and by the Joint Services Electronics Program. 3
1047-3203/90 $3.00 Copyright 0 1990 by Academic Press, Inc. All rights of reproduction in any form reserved.
4
ROGER W. BROCKETT
and the vector tion of a plane. curve lying on the underlying not.
fields associated with perspective projecThey show that in this context a biquartic a plane has enough variation in it to reveal flow but that curves of lower degree do
2. PRELIMINARIES
ON LEAST SQUARES
The mathematical idea which is basic in our search for robustness, one which is standard in many areas, is nothing more than a suitable version of the least-squares method of curve fitting. Suppose that we have a set of observations depending on a parameter s:
I: II Y,(S)’ Y2b)
Y(S) =
*
Y,(S).
Suppose further that these observations are generated by a linear transformation acting on an unknown, but constant, vector
x zzs-1 P TT(s)y(s)ds. ia If S is singular then the minimizing value of x is not unique. The invertibility of the symmetric, nonnegative definite matrix S is a necessary condition for this problem to have a unique answer. (See [14, p. 711.) From a geometrical point of view, least-squares problems are solved by constructing a suitable orthogonal projection. The symmetric matrix S is diagonalized by an orthogonal transformation which sets up the orthogonal projection; i.e., it orthonormalizes the columns of T with respect to the inner product structure defined by the quadratic form in error criterion. Such matrices arise in any systematic approach to orthonormalization (e.g., the Gram-Schmitt process) and are often called gramians. Note that if the weights on the various components of the error are changed then the gramian is pre and postmultiplied by a suitable digonal matrix: S -+ DSD. Increasing the level of generality slightly now, suppose we want to minimize e = I : (T(t)x - y(t), T(t)x - y(t))& subject to constraints
Xl
on x of the form Px = 6.
x2
x=
.
Xfl
y(s)
= T(s)x,
a! 5 s 5 p.
Since this equality can be thought of as giving one equation for each s it is generally an overdetermined system and we cannot expect to solve it. However, if we seek the x which minimizes the square of the difference between the left- and right-hand sides, i.e.,
To avoid unnecessary complications, we assume that the constraints are self-consistent and expressed in a minimal way-i.e., that P is p by n and of rank p. In this case we see that for S = I” TT(t)T(t)dt a we have x = s-1
I a’ TWy(t)dt
with A being given by the solution of PS-’ I”01TT(t)y(t)dt
then there is a well-known the following lemma.
formula for x which we give in
LEMMA 1. Zf the matrix
s = I,” TT(s)T(s)ds is nonsingular
then the value of x which minimizes 7 is
+ PTA,
+ PPTh = b,
i.e., by A = (PPT)-‘(b
- PS-’ I,” TT(t)y(t)dt)
.
3. VECTOR FIELDS AND GRAMIANS
In this section we formulate which involves
and answer a question
GRAMIANS,
INVERSES,
(i) A finite collection of vector fields {Fi}y=r which have a common domain of definition in euclidean space. (ii) A distinguished subset S of the domain of definition of the Fi. (iii) A “directional” vector field T defined on S. Remark. There are two prototype situations which have been considered in the literature and which are part of the motivation for our setup. In one case S is the image of a surface with reflectance R(x, y) and T is (aR/ax, aRl8y). In the other case S is the image of a contour and T is the unit normal to the contour. In either case the set of vector fields would include the infinitesimal generators of translation in the horizontal and vertical, rotation, and possibly other vector fields as dictated by the image formation law and the surface shape. Problem. Suppose that on S we are given the dot product between T and an unknown vector field FO. Find al,a2,. . . , a,, such that the linear combination a I F1 + a2F2 + ... + a,F,, approximates F0 in the sense that
q = Is((FO- 2 aiF;,T))2dS is as small as possible. Here (F, T) denotes the ordinary dot product of the vectors F and T. The subset S will either be the union of the closure of a finite collection of open sets where dS indicates the ordinary volume measure or else S will be the union of a finite set of smooth curves where dS indicates the differential of arc length. As can be seen from the previous section, the basic ingredient in the solution of this least-squares problem is the gramian G whose 0th entry is gij = I s (Fi, T) * (4,
T)ds.
Applying the results of the previous section we can assert the following. 2. Zf the gramian G is nonsingular, then the optimal value of a = (at, a2, . . . , a,,) is unique and is given in terms of (FO, T) by LEMMA
a = G-t I s
VI, 0 [1 . . .
OPTICAL
5
FLOW
fields in the plane since this is the situation which arises in computer vision. 4.
THE VECTOR
FIELDS
OF OPTICAL
FLOW
As an object moves through space, its image flows across the image plane. If (x, y) is used to denote the horizontal and vertical coordinates, respectively, of a point on the image plane, then the velocity vector associated with the image, expressed as a function of x and y, defines a vector field on part of the image plane. There are a number of formulations of the image formation process including orthographic and projective. Rigid motions may be general three-dimensional motion or may be constrained in some way. The surfaces that rigid objects present to the camera may be planar, quadratic, or of some other form. Each combination of possibilities gives rise to vector fields on the image plane having a particular form. We illustrate using some well-known examples. 4.1. Rigid 20 Motion, Motion Plane
Zmage Plane Parallel to
This is a very simple situation. We consider a planar rigid object moving in a plane parallel to the image plane. Because the plane of motion is taken to be parallel to the image plane, all points in the motion plane are at equal depth and it is immaterial whether orthographic or projective image formation is assumed. Rigid motion in the plane is characterized by two translational velocities and one rotational velocity. In this case the velocity in the image plane is given by making a choice of CX,/3, w in the equation.
[r:J=I”
rl[:l+[;].
4.2. Rigid 30 Motion, Planar Surface
Orthographic
Projection,
In this case, it is convenient to choose a coordinate system for three-dimensional euclidean space with one axis, which we take to be the Z axis, perpendicular to the image plane. Orthogonal projection is represented by POW, Y, Z) = (X9 Y) = (x,
Yh
Wo, TW.
(Fn, T>
Zf G is not invertible, unique.
AND
then the optimal
value of a is not
Note that for the purposes of the lemma the dimensionality of the euclidean space is irrelevant; however, in the remainder of this paper we use only the case of vector
As a rigid body moves in three-dimensional space, a point vi fixed in it transforms by a formula of the type
@vi + b = sl) with 0 being an orthogonal matrix. Associated with a time-dependent euclidean transformation are a pair of three vectors, the rectilinear velocity
6
ROGER
W. BROCKETT
and the angular velocity. The angular velocity is intrinsic whereas the rectilinear velocity must be defined with respect to a choice of origin. As is well known, associated with a time-dependent orthogonal matrix o(t) is an angular velocity matrix
[ I 0
03
R=
dependent with
02
-010
on the angular velocity
0
vector (0, , 02, 03)
$ o(t) = n(t)@(t). Affine transformations form
$13
023
Y’
.
-w2
01
-w3
+[Ihbz+[I z"=1TX'11 Note that the problem of determining the full three-dimensional motion from its two-dimensional flow cannot be separated from the identification of the quantities p, q, and Z. which relate to the shape and depth. The key point here is that by differentiation of this formula one sees that the rigid motion of a plane in three-dimensional space gives rise to an affine velocity field on the image plane. To establish this fact, note that in terms of CR, as introduced above, we have
are more general. They take the
with 9 being any invertible matrix. In two dimensions there :s a three-parameter family of euclidean transformations and a six-parameter family of affine transformations. The corresponding numbers in three dimensions are 6 and 12. Orthogonal projection ignores the 2 axis, and thus the distance between points on a rigid body as projected onto the image plane will generally change as the body undergoes a rigid three-dimensional transformation. In fact, we see that POqf and POqi are related by
If we evaluate this for 0 equaling the identity matrix we get a differential version of the projection equation which involves the rotational and translational velocities
[;I=[-13 :][;I+[z:J [I +[;;]+[;]zo. X
PoOq; + POb = P&,
and thus while the coordinates of the projected image after transformation are linearly dependent on those before transformation, they also depend on the unknown Z component. If (X, Y, Z) denotes a point on the surface of an opaque object, a formula for Z in terms of X and Y is tantamount to knowing the shape of the object. Thus the above expression requires a knowledge of the shape in order to be empirically useful. If the visible surface of the object is a plane, it is described by an equation Z = pX + 4 Y + ZO. In this case rigid motion in three-dimensional space yields a (possibly singular) affine motion in the image plane. Writing
we can express this transformation
as
LPY
41
Y
The right-hand side is, indeed, an affine function of x and y. It is to be noted that the set of all possible affine functions of two variables is a six-parameter set. 4.3. Rigid 30 Motion, Planar Surface
Perspective
Projection,
As a third example of planar vector fields generated by optical flow, we consider the optical flow generated by rigid 3D motion of a planar surface as above but now with a perspective image formation law. That is to say, a point (X, Y, Z) in E3 is projected onto the image plane according to P(X, Y, Z) = (X/Z,
Y/Z) = (x, y).
GRAMIANS,
Differentiation
INVERSES,
yields v, = klZ VY
- XilZ2
= ivz - YilZ2.
Assuming that the visible surface is a plane, say Z = pX + qY + ZO, it has been shown that on the image plane v, and v?, are given by quadratic functions of x and y: v, = al + a2x + u3y + @X2 + ujxy + a6y2
AND
OPTICAL
7
FLOW
where s is arc length and IZis the unit normal to the curve. There is no need for the boundary curve to be connected. The integral may include contributions from many arcs. We propose to call this the shape grumian. Feature point tracking is not included explicitly in this formalism but can be thought of in the following way. If an object has feature points p1 , p2, . . . , pk which can be tracked, surround each of them with a small circle having radius r and center pi. Evaluate the shape gramian associated with these k circles. The limit
v, = b, + b2x + b,y + b4x2 + b5.q + bby*.
The interested reader can consult Waxman and Wohn [6] for a derivation. (Note that there are but eight free parameters.) 5.
THE DIRECTIONAL
VECTOR
FIELD
Horn and Schunck [2] and Wohn et al. [15] have reported results on using the “total” derivative of the intensity f
Z(x, y) = R,i
+ Ryj
as a means for determining ,i and j, (see, e.g., [6]). This makes sense if the object in question has a significant change in its intrinsic reflectance and if the ambient lighting is sufficiently uniform so that the convected invariance assumption is appropriate. To fit this into our setup we take S to be the image of the object in question and we take T to be the vector field T = VI
(V = gradient).
In this case the formula for the relevant gramian is go = I, (Fi, V R)(Fj,
V R)dxdy.
We propose to call this the reflectance grumiun. In the work of Marr and Ullman 111, Hildreth [3], and Waxman and Wohn [6] the problem of determining the motion of an opaque object based on the movement of its boundary is discussed. It is argued that in such a situation it is the component of the velocity which is normal to the boundary which can be determined. In this case the set S is the boundary of the object in question and T is the unit normal to the boundary. This is illustrated in Section 7. The expression for the gramian is
gu = I s (F;, n)(Fj, n)ds,
is finite and depends only on the location of the feature points. We propose to call this thefeature point grumiun. Of course it may happen that some feature points are more reliable than others. If this can be quantified, then the above definition can be modified by associating a different scaling of the radii, depending on the robustness of the associated point. These three cases are to be thought of as prototypes. There will be situations in which both reflectance and shape information are significant and situations in which feature point information needs to be combined with shape and reflectance information. In these cases the appropriate gramian may contain contributions from V Z over areas, the outward pointing normal vector over arcs, and point contributions as well. One principal motivation for formulating the problem of optical flow recovery in this generality is to be able to treat such hybrid problems. 6.
THE AFFINE
REFLECTANCE
GRAMIAN
The general principles discussed above will now be applied to some specific settings. Obviously not all possibilities can be worked out here, so we have chosen to look in some detail at a few standard problems. The problem of determining the best affme approximation to optical flow based on shape data is treated in this section. Let I, be a one-parameter family of simple closed curves in the plane. Assume that these curves are represented as a function of t (= time) and their arc length s according to
x0, 4 = $0, s) YG,
s)
=
40,
s),
with $ and 4 being piecewise differentiable functions oft and s. Let R(x, y) denote the intrinsic reflectance of an object which projects onto the image plane and boundary Ir. Let Z be the perceived intensity on the image plane.
8
ROGER
W. BROCKETT
tor fieldfr(x, y) = (YX + by + e, &4x, y) = cx + dy + J Explicitly, we assume that Z is given and we want to find a, b, c, d, e, f such that
Let F be a smooth vector field on the image plane. Consider the expression for the derivative along F $ Z(x, y) = R,Ec + Ryj
r)=
For i = fi(x, y); j, = &(x, y) we have
i=R,fi
ARG=I\
is minimized. Define a 6 x 6 symmetric matrix AZ?G by the expression (the integration is over the interior of F)
R2x2 x
R:XY
R,Ryx2
R,R,xy
R:x
R,R,x
R;xy
R:y2
R,R,xy
R,R,y2
R:x
RxRyx
R,Ryx2
R,R,xy
R;XY
R?
R,R,x
R,R,x
RxR,xy
R,R,y=
R;XY
R;
RxR,y
RyRyy
R;x
R:x
RxRyx
R,R,x
R,R,x
R;x
THEOREM 1. Zf ARG is nonsingular, (a, 6, c, d, e, f) which minimizes 11 is
then the value of
zxx LY ZYX
ZYY 1.X
the minimizing
R,R,
R::
R,Ry
Rf
dxdy.
lem if ARG is nonsingular. The nonsingularity of ARG and, in fact, the size of the smallest eigenvalue of ARG, influences the robustness of this solution. We will call ARG the afine rejiectance gramian associated with Z and R. Later on we will discuss in qualitative terms the nonsingularity of ARG.
7. THE AFPINE SHAPE GRAMIAN Suppose that Z(x, y) is a function which is one on R and zero outside R. Such a case corresponds to a silhouette moving across a plain background. All the information on the motion comes from the moving edge. Let (n,. . ,., n,),, be the components of the outward pointing normal to F at s. This leads to a formula involving the symmetric matrix
value of (a, b, c, d,
This result is an immediate consequence of the results cited above. It provides a complete solution of the prob-
1
n2x2 Y
AGS = 1
R:
R&Y
-
k, y)drdy.
ZY Zf ARG is singular, e, f) is not unique.
+ by + e)
+Ryf2.
The first problem we consider is that of associating to of F by an affine vec-
11,
R (i - R&x
- R,(cx + dy + f))2 dxdy
R the least-squares approximation
= (ARG)-’
II
+y
nxnyx
2
nxnyw 2
4xy
n?y2
nxnyxy
wyx2
nay-y
n2x2
d-v
nxnyxY
nxnyY
nZxy
n2y2
nix
nix
nay
vy
n2w
+
wyx
wyx
vyY
nfx
nxnyx
n,2x
w,x
nay
+
w,y
4
vy
ds
4y 2
nY
:
wy 4
2
I
GRAMIANS,
INVERSES,
which we call the u~@e shape gramian. The leastsquares fit thus results in a formula for the motion parameters . THEOREM 2. Zf the ASG is nonsingular then the value of (a, b, c, d, e, f) which minimizes 7) is
AND
OPTICAL
9
FLOW
From the above congruent formula, we see that the shape gramian of a parallelogram is also nonsingular since affine transformation will take a square into a parallelogram. On the other hand, the ASG of an ellipse is singular since that of a circle is singular and an affine transformation can be found to convert an ellipse into a circle.
w nYX
= (ASG)-’ 1,.
w
8.
v,(s)ds.
TWO-DIMENSIONAL RIGID MOTION THE IMAGE PLANE
PARALLEL
TO
w nY n, I-
Zf ASG is singular the minimization not unique.
of (a, b, c, d, e, f) is
This formula is the appropriate one when only the edges or silhouette of an object can be observed and no shading information is available. Since the ASG can be difficult to evaluate, it is useful to discuss how it depends on coordinate transformations. The key point is that it transforms by a congruence which depends on the underlying transformation of coordinates in a simple way. THEOREM 3. Zf ASG is the ufJine shape grumiun with respect to a coordinate system (x, y) and if (x’, y’) is related to (x, y) by
[t:]=[:iI [:1+[:19 then the ufJine shape grumian with respect to (x’, y’) is given by PTSP, with P being nonsingular if and only if the mapping from (x, y) to (x’, y”) is invertible. It is possible to express P quite briefly in terms oftensoral transformation laws which make the naturality of this construction more explicit. In order to fix ideas we evaluate the affine shape gramian for a square, chasing a coordinate system centered at the center of the square and aligned with the edges. If the square is two units on a side the ASG is
This section is a digression devoted to the study of the special case of two-dimensional euclidean motion. This is of interest in its own right and as a simplification of the affine case treated above. Consider a curve in two-dimensional euclidean space which is given by {x, y[ f (x, y) = O}. Pick an orientation for the curve and let s denote arc length along the curve. Because the group of rigid motions in the plane is threedimensional and no projection is involved, we associate with such a curve a 3 x 3 symmetric matrix which is the appropriate one for extracting the translational and rotational velocities associated with rigid motion. The process involves two steps. First, we introduce a construction which involves an arbitrary choice origin. In step two we will eliminate this dependence on how the curve is displaced from the origin by associating to the curve a standard choice of origin. If the curve undergoes a rigid motion caused by an xvelocity ,, a y-velocity y, and a rotation 8 about the origin, then at a point (x, y) on the curve, the curve moves in its normal direction with a velocity
l L”]*([:]+[1,lq vw=l@p fy That is, the normal velocity is 1
EXAMPLE.
400000 040000
1
000040 00000Y
v(s) = *
+ f;
.((fc 3f, 3 Yfx
- xfyh t-k j, f%>.
This equation is of the form 4s) = (P(S), i).
1
1
If we regard the various values of v(s) as given then we can ask how best to recover the values of the parameters ,, y, 6 which characterize the rigid motion. Assuming we have scaled x and y so that it is reasonable to minimize the sum of the squares of the errors in ,, y, and 0, then the results of Section 2 apply. In this case we have
10
ROGER
W. BROCKETT
not determine i, j, 8, uniquely. It is not difficult to see that S(F) is singular if and only if I’ has constant curvature. That is, I must be a line segment or an arc of a circle. The choice of coordinate frame for this calculation has, until now, been arbitrary. We should, thus, ask how S(F) changes if the coordinate system is changed? Say we move the origin and rotate coordinates. Introduce
with
[1 PI(S)
ESG
=
Ir
p2W
[PI@),
p2bh
mWld~.
P3b)
Y’
We call ESG the two-dimensional euclidean shape gramian of the curve relative to the given frame F. If it is singular then the normal component of the velocity does
cos I30 S(F) =
If we partition
::sHBu] [:I+[I. [ I=[::iYio X’
Going back to the definition we see that this gives a new s:
-sin 80 cos 130 y0
y.
S as
rd
s = I”:,,1
a
Ls then it becomes
S=
1
=r @TOT
L(OTa + Om)
(OTa + Om)T r - 2md + aTTa
Thus the determinant and trace of Tare intrinsic. Moreover, since there is always an a such that Ta + m = 0 we can always choose an origin such that the ESG is diagonal. It is not too difficult to see that except for curves whose normal never changes this choice is unique. For straight line segments this constrains the origin to the locus of points defined by the perpendicular bisector of the segment. We refer to this choice of origin as the centroid of the curve and refer to the axes which diagonalize T as the symmetry axes of the curve. EXAMPLE 1. Consider a straight line of length r. If the origin of the coordinate system is taken to be the midpoint of the line and if the line has slope CYthenfr(s) = xl-, f2(s) = ax/V?%?, and s goes from -r/2 to r/2. In this case
1
x0 T
sin 80 ~0’
.
a
s
1
-03
-as
s2
1 ds
1 a2 ff 1 s=r[(y 1 [
a2/a2 + 1 arta + 1
0
(r/a2 + 1
l/O? + 1
0
0
0
r2/12
.
Note that in this case S is singular. If we move the origin to the endpoint of the line then S becomes
a2 + 1
1
1
r/2
r/2
r2/3
,
which is, of course, also singular. EXAMPLE 2. Consider an equilateral triangle with coordinate system at the centroid of the triangle. In terms of the perimeter p, S is given by
LO
0
p/271
GRAMIANS,
INVERSES,
EXAMPLE 3. A regular n-gon of perimeter p has a centroid located at its center and has gramian shape
s=
-p/2
0
0
0
p/2
0
-Q
0
p3/3n2
I-
We note that the ESG of a plane curve expressed with respect to the natural coordinate system is invariant with respect to rigid motion and thus can be expressed in terms of the arc length and the curvature, alone. For example, the trace of T, defined above, is the arc length.
9. CONCLUSIONS
Relative motion between the objects being viewed and the camera yield a changing pattern on the image plane. Depending on the situation, the piecewise smooth flow can be represented as a member of a parametrized family of vector fields of a specific type. Solving the inverse problem, i.e., the problem of finding the relative motion in three-dimensional space given the changing patterns in the image plane, often involves the computation of a representation of the two-dimensional flow field in the image plane. We have shown here that there is a natural common framework for this process. It involves a generalized inverse of an integral operator and through the mathematical formulation of this process we are lead to necessary and sufficient conditions for the determination of the optical flow.
AND OPTICAL
II
FLOW
REFERENCES 1. D. Marr and S. Ullman, Directional selectivity and its use in early visual processing, Proc. Roy. Sot. London Ser. B. 211, 1981, 151180. 2. B. K. Horn and B. G. Schunck, Determining optical flow, ArtiJicial Intelligence 17, 1981, 185-208. 3. E. C. Hildreth, Computations underlying the measurement of visual motion, Art$cial Intelligence 23, 1984, 309-354. 4. R. J. Schalkoff, Image Motion Analysis Using the Concept of Weak Solutions to Distributed Parameter Systems, CVPR, 1983. 5. K. Kanatani, Structure from motion without correspondence: General principle, Proc. IJCHA, 1985, 886-888. 6. A. M. Waxman and K. Wohn, Contour evaluation, neighborhood deformation and global image flow: Planar surfaces in motion, Internat. J. Robotics Res. 4, 1985, 95-108. 7. J. J. Koenderink and A. J. Van Doom, Local structure of movement parallax of the plane, J. Opt. Sot. Amer. 66, 1976, 717-723. 8. H. C. Longuet-Higgins and K. Prazdny, The interpretation of a moving retinal image, Proc. Roy. Sot. London Ser. B. 208, 1980, 385-397. 9. R. Y. Tsai and T. S. Huang, Uniqueness and estimation of threedimensional motion parameters of rigid objects with curved surfaces, in Proceedings, PRIP 1982, Las Vegas, NV, pp. 112-118. IO. S. Ullman, The Interpretation of Visual Motion, MIT Press, Cambridge, MA, 1979. Il. A. M. Waxman and S. Ullman, Surface structure and 3D motion from image flow: A kinematic analysis, Znternnt. J. Robotics Res. 4, 1985, 72-94. 12. A. R. Bruss and B. K. P. Horn, Passive navigation, Comput. Vision, Graphics Image Process. 2, 1983, 3-20. 13. A. Bandypopadhyay and J. Aloimonos, Perception of rigid motion from spatio-temporal derivatives of optical flow, Unpublished manuscript . 14. R. W. Brockett, Finite Dimensional Linear Systems, Wiley, New York, 1970. 15. K. Wohn, L. S. Davis, and P. Thrift, Motion estimation Based on multiple local constraints and nonlinear smoothing, Pattern Recognition 16, 6, 1983, 563-570. 16. E. J. DeIp and P. H. Eichel, Motion parameter estimation: Developing the scale change parameter, Unpublished manuscript.