Method of modifying a time-varying image sequence by estimation of velocity vector fields

Abstract

A method of modifying a time-varying image sequence in which a parametric velocity vector field that characterizes changes in successive images of the image sequence is estimated. The estimating includes combining optical flow constraints, directional smoothness constraints and regularization constraints in a functional of the estimate of the parametric velocity vector field such that the weighting coefficients of the functional are functions of the parametric velocity vector field to be estimated. A system of nonlinear equations that arise from the optimality criterion of the functional are then solved. Once the system of nonlinear equations is solved, the estimate of the parametric velocity vector field is applied to modify at least one image in the image sequence.

Description

FIELD OF THE INVENTION

The present invention relates to the field of image processing, and more particularly, to a method of modifying a time-varying image sequence based on an estimate of the velocity vector field.

BACKGROUND OF THE INVENTION

A time-varying image sequence is a sequence of images of a given scene with each successive image taken some time interval apart from the one preceding it. The velocity vector field of an image is defined as the collection of the two-dimensional velocity vectors which are the projections of the three-dimensional velocity vectors of the visible points in the scene. If the scene being imaged changes gradually with time, and if the changes are mostly due to the relative movements of the physical objects in space, then the corresponding changes in the successive images of the sequence can be used to estimate the velocity vector fields of the images.

Reliable estimation of velocity vector fields is very important for the analysis of time-varying image sequences. The temporal variations in the images of the sequence specified by the velocity vector fields can be used to extract spatial information about the scene. They can also be used to extract information characterizing the movements of the objects in the scene. Stereopsis can be considered as a special case of image-sequence analysis where the number of images is restricted to two and the displacement direction is known.

There are two principal approaches to the problem of estimation of the velocity vector field: the feature-based matching approach and the spatio-temporal gradient approach. For feature-based matching approach, the image points with significant variations in the values of the time-varying image function, called feature points, are identified in both images. The estimation of the velocity vector field is accomplished by matching the feature points of one image to the feature points of the other image. The spatio-temporal gradient approach is based on the constraint imposed on each velocity vector relating the spatial gradient of time-varying image function to the temporal derivative of the time-varying image function.

The spatial variation in the time-varying image function utilized in the above approaches do not provide sufficient information to determine the estimate of the velocity vector field. In the feature-based matching approach the velocity vectors can be estimated only on a sparse set of image points, while in the spatio-temporal gradient approach at most one constraint is imposed on two components of each velocity vector. To overcome these difficulties, it has been proposed, Horn and Schunk, "Determining Optical Flow", Artif. Intell. 17, 1981, pp. 185-203, that velocity vector fields should vary smoothly from point to point on the image plane. This requirement enabled estimation of both components of the velocity vector at each image point; however, it forced the estimate of the velocity vector field to vary smoothly across the occluding boundaries. Several approaches, which are based on the selective application of the smoothness requirement, have been proposed to overcome this difficulty.

Yachida used estimates of the velocity vectors at prominent feature points as reliable initial estimates and sequentially propagated them into neighboring image points. Davis et al first computed estimates of the velocity vectors at corner points and then propagated them along the contour lines between corner points. Hildreth estimated the velocity vector field along zero-crossing contours. Nagel used image points with significant second-order spatial variation to estimate of velocity vector fields. Nagel and Enkelmann investigated oriented smoothness constraints for the estimation of velocity vector fields. Terzopoulos introduced controlled continuity measures to overcome the discontinuities in visual problems. Hierarchical approaches to the problem of estimation of velocity vector fields have also been investigated.

There is a need for a method of modifying a time-varying image sequence based on an improved method of estimating the velocity vector field.

SUMMARY OF THE INVENTION

This and other needs are met by the present invention which provides a method of modifying a time-varying image sequence in which a parametric velocity vector field that characterizes changes in successive images of the image sequence is estimated. The estimating includes combining optical flow constraints, directional smoothness constraints and regularization constraints in a functional of the estimate of the parametric velocity vector field such that the weighting coefficients of the functional are functions of the parametric velocity vector field to be estimated. A system of nonlinear equations that arise from the optimality criterion of the functional are then solved. Once the system of nonlinear equations is solved, the estimate of the parametric velocity vector field is applied to modify at least one image in the image sequence.

In another embodiment of the present invention, the estimation of the parametric velocity vector field includes combining optical flow constraints, directional smoothness constraints and regularization constraints in a functional of the estimate of the parametric velocity vector field such that weighting coefficients of the functional are functions of the parametric velocity vector field to be estimated; and providing a dual feed-back interaction between the optical flow constraints and the directional smoothness constraints.

These and other objects of the present invention will become more apparent when taken in conjunction with the following description and drawings wherein like characters indicate like parts and which drawings form a part of the present description.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 illustrates a schematic block diagram of apparatus useable to perform the method of the present invention.

FIG. 2 illustrates a flow chart for a method of modifying an image in a time-varying image sequence in accordance with an embodiment of the present invention.

FIG. 3 illustrates a flow chart for a method of modifying an image in a time-varying image sequence in accordance with another embodiment of the present invention.

FIG. 4 illustrates a flow chart for a method of modifying an image in a time-varying image sequence in accordance with another embodiment of the present invention.

FIG. 5 illustrates a flow chart for a method of modifying an image in a time-varying image sequence in accordance with another embodiment of the present invention.

FIG. 6 illustrates a flow chart for a method of modifying an image in a time-varying image sequence in accordance with another embodiment of the present invention.

FIG. 7 illustrates a flow chart for a method of modifying an image in a time-varying image sequence in accordance with another embodiment of the present invention.

FIG. 8 illustrates a flow chart for a method of modifying an image in a time-varying image sequence in accordance with another embodiment of the present invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

In implementing the method of the present invention, conventional physical image processing equipment is used. This equipment would include, as seen in FIG. 1, for example, a camera 10, an image digitizer 12, a computer 14, and a display 16 that are entirely conventional.

The following is a description of the formulation of the estimation problem. The images are of a three-dimensional scene and their spatial partial derivatives are take within some time interval T. The initial images of a three-dimensional scene are, in general, discontinuous functions and their spatial partial derivatives are not defined at the points of discontinuity. The points of discontinuity are, often, among the most important points in the images and cannot be easily ignored. To overcome these difficulties the initial images are treated as generalized functions and their partial derivatives as generalized partial derivatives. These generalized functions are defined on some subset of the set of infinitely differentiable functions. The functions from the subset are called "testing functions". Parametric families of secondary images are then introduced through evaluations of the generalized functions associated with the initial images on the specific parametric family of functions taken from the set of testing functions. The secondary images are infinitely differentiable, and their partial derivatives can be obtained through evaluations of the generalized functions associated with the initial images on the partial derivatives of this specific parametric family of functions. This process can be described as follows.

One way of constructing an initial image if by projecting the light reflected by the objects in the scene at a given moment of time t.epsilon.R (here R is a one-dimensional Euclidean space) onto a two-dimensional Euclidean space R.sup.2, called the "image plane", and then identifying the irradiance value .xi.(x,y,t) of each point (x,y) in the image plane R.sup.2. The function .xi.(x,y,t) defined as above will be called the "irradiance image function". The value .xi.(x,y,t) of the irradiance image function at the point (x,y).epsilon.R.sup.2 and the time t.epsilon.R is assumed to be roughly proportional to the radiance of the point in the scene being imaged, which projects to such a point (x,y) at the time t for every (x,y).epsilon.R.sup.2 and every t.epsilon.R. Different irradiance image functions .xi.(x,y,t), (x,y,t).epsilon.R.sup.3 (here R.sup.3 is a three-dimensional Euclidean space) can be obtained by changing the following aspects of the image formation process: the direction of a light source illuminating the scene, the color of such a light source, and the spectral responsivity function that is used to compute the irradiance.

Assume that each irradiance image function .xi.(x,y,t) is locally integrable with respect to the Lebesgue measure dx dy dt in R.sup.3 and thereby can be used to form a continuous linear functional .GAMMA..sub..xi. (generalized function) defined on the locally convex linear topological space .PHI.(R.sup.3). The space .PHI.(R.sup.3) consists of all infinitely differentiable functions having compact supports in the set R.sup.3. This means that for each function .phi..epsilon..PHI.(R.sup.3) there exists a closed bounded subset S.sub..phi. .epsilon.R.sup.3 such that the function .phi. is equal to zero at the points that are outside the subset S.sub..phi.. The topology of the space .PHI.(R.sup.3) is defined by a certain family of semi-norms. The functions .phi. from the set .PHI.(R.sup.3) will be called the "testing functions". The value of the generalized function .GAMMA..sub..xi. associated with the irradiance image function .xi. at the testing function .phi..epsilon..PHI.(R.sup.3) is defined by the following relation: ##EQU1## The generalized function .GAMMA..sub..xi. associated with different irradiance image functions .xi.(x,y,t) are united into the family {.GAMMA..sub..xi. .vertline..xi..epsilon..XI.}. Here, to simplify the notation, the symbol .xi. takes the role of the index (parameter) specifying the generalized function .GAMMA..sub..xi. associated with a particular irradiance image function .xi.(x,y,t) in addition to its role as a part of the notation for such as irradiance image function. The symbol .XI. denotes the set of indices.

Another way of constructing an initial image at a given moment of time t.epsilon.R is by identifying the set of points M(t).OR right.R.sup.2 in the image plane, called the "feature points", where significant variations in the projected light patterns represented by the irradiance image functions takes place at the time t, and then assigning feature value .eta.(x.sub..mu.,y.sub..mu.,t.sub..mu.) to each feature point (x.sub..mu.,y.sub..mu.,t.sub..mu.).epsilon.M(t). The set of feature points M.ident..orgate./t.epsilon.R M(t) is assumed to be a closed subset of the space R.sup.3. The function .eta.(x.sub..mu.,y.sub..mu.,t.sub..mu.) defined on the set M as above will be called the "feature image function". Different feature image functions .eta.(x.sub..mu.,y.sub..mu.,t.sub..mu.), (x.sub..mu.,y.sub..mu., t.sub..mu.).epsilon.M can be obtained by changing the criteria for selecting the set of feature points M and the criteria for assigning the feature values .eta.(x.sub..mu.,y.sub..mu.,t.sub..mu.). The set M may, for example, be a finite combination of the following four types of the subsets of the space R.sup.3 : the three-dimensional regions, the two-dimensional regions, the one-dimensional contours, and the isolated points.

By the family B.ident.B(M) of "Borel subsets" of the set M shall be meant the members of the smallest .sigma.-additive family of subsets of M which contains every compact set of M. Let .mu.(B) be a .sigma.-finite, .sigma.-additive, and real-valued measure defined on the family B of Borel subsets of the set M. The feature image function .eta.(x.sub..mu.,y.sub..mu.,t.sub..mu.), (x.sub..mu.,y.sub..mu.,t.sub..mu.).epsilon.M is assumed to be .mu.-measurable measurable and thereby can be used to form a continuous linear functional .GAMMA..sub..eta. (generalized function) defined on the locally convex linear topological space .PHI.(R.sup.3). The value of the generalized function .GAMMA..sub..eta. associated with the feature image function .eta. at the testing function .phi..epsilon..PHI.(R.sup.3) is given by the following relation ##EQU2## The generalized functions .GAMMA..sub..eta. associated with different feature image functions .eta.(x.sub..mu.,y.sub..mu.,t.sub..mu.) are united into the family {.GAMMA..sub..eta. .vertline..eta..epsilon.H}. Here, to simplify the notation, the symbol .eta. takes the role of the index (parameter) specifying the generalized function .GAMMA..sub..eta. associated with a particular feature image function .eta.(x.sub..mu.,y.sub..mu.,t.sub..mu.) in addition to its role as a part of the notation for such a feature image function. The symbol H denotes the set of indices.

Given a generalized function F defined on the locally convex linear topological space .PHI.(R.sup.3), and given non-negative integer constants m.sub.x, m.sub.y, m.sub.i the "generalized partial derivative" ##EQU3## of the generalized function F is the generalized function ##EQU4## defined on the locally convex linear topological space .PHI.(R.sup.3) as follows. The value of the generalized function ##EQU5## at a testing function .phi..epsilon..PHI.(R.sup.3) is the value of the generalized function F itself at the testing function ##EQU6##

A "combined generalized initial image function" will be defined as a linear combination of the generalized partial derivatives of the generalized functions associated with the irradiance image functions and of the generalized partial derivatives of the generalized functions associated with the feature image functions. Let .lambda..ident.{.lambda..sub..xi., .lambda..sub..eta. .vertline..xi..epsilon..XI., .eta..epsilon.H} be a set of real-valued constants; let g(.lambda.) be an index attached to the set .lambda.; and let m.sub.g(.lambda.) .ident.{m.sub..xi.,x, m.sub..xi.,y, m.sub..xi.,t, m.sub..eta.,x, m.sub..eta.,y, m.sub..eta.,t .vertline..xi..epsilon..XI., .eta..epsilon.H} be a set of non-negative integer constants corresponding to the set of constants .lambda.. Then the combined generalized initial image function corresponding to the set of constants .lambda. is the generalized function .GAMMA..sub.g(.lambda.) defined on the locally convex linear topological space .PHI.(R.sup.3). The value of the combined generalized initial image function .GAMMA..sub.g(.lambda.) corresponding to the set of constant .lambda. at the testing function .phi..epsilon..PHI.(R.sup.3) is given by the relation ##EQU7## and is called the "observation" of the combined generalized initial image function .GAMMA..sub.g(.lambda.) corresponding to the set of constants .lambda. on the testing function .phi.. The combined generalized initial image functions .GAMMA..sub.g(.lambda.) corresponding to sets of constant .lambda. with different values are united into the family {.GAMMA..sub.g(.lambda.) .vertline..lambda..epsilon..LAMBDA.}. Here the symbol .LAMBDA. denotes the family of all different sets of constant .lambda., while the symbol g denotes the one-to-one mapping from the family .LAMBDA. onto the set of indices denoted by the symbol G. In what follows, to simplify the notation, the argument .lambda. appearing in the notation for a combined generalized initial image function is omitted and the symbol .GAMMA..sub.g, g.epsilon.G is used instead of the symbol .GAMMA..sub.G(.lambda.), .lambda..epsilon..LAMBDA. to denote it.

Assume that there is a given fixed testing function .PSI..epsilon.(.PHI.(R.sup.3), which will be called the "measurement function". (An example of the measurement function .PSI. is given later.) Let g.epsilon.G; then for every point (x,y) of some convex bounded subset .OMEGA..OR right.R.sup.2 (such points will be called "image points"), for every moment of time t.epsilon.T, and for every value of the parameter .sigma..epsilon.[1, .varies.), the value of the component g.sup..sigma. (x,y,t) of the image corresponding to the parameter value .sigma. at the image point (x,y) and time moment t is determined as the observation .GAMMA..sub.g (.PSI..sub.x,y,t.sup..sigma.) of the combined generalized initial image function .GAMMA..sub.g on the testing function .PSI..sub.x,y,t.sup..sigma. .epsilon..PHI.(R.sup.3), ##EQU8## Note that, to simplify the notation, the symbol g is used as the part of the notation g.sup..sigma. for the component of the image corresponding to the parameter value .sigma. in addition to its role as the index for such a component.

The parametric vector valued function g.sup..sigma. (x,y,t).ident.{g.sup..sigma. (x,y,t).vertline.g.epsilon.G}, (x,y,t).epsilon..PSI.xT, .sigma..epsilon.[1, .varies.) with each component g.sup..sigma. (x,y,t), g.epsilon.G defined on the set .PSI.xT as the observation .GAMMA..sub.g (.PSI..sub.x,y,t.sup..sigma.) of the combined generalized initial image function .GAMMA..sub.g on the testing function .PSI..sub.x,y,t.sup..sigma. .epsilon..PHI.(R.sup.3) specified by the relation (2-4) for every parameter value .sigma. from the set [1, .varies.) (here the symbol .OMEGA.xT denotes the Cartesian product of the set .OMEGA. on the set T) will be called the "parametric time-varying image function". The points (x,y,t) from the set .OMEGA.xT will be called "time-varying image points". By the "image" corresponding to the parameter value .sigma..epsilon.[1, .varies.) and taken at the time moment t.epsilon.T is meant the collection, {g.sup..sigma. (x,y,t).vertline.(x,y).epsilon..OMEGA.}, of the values of the parametric time-varying image function g.sup..sigma. (x,y,t), (x,y,t).epsilon..OMEGA.xT. The collection {g.sup..sigma. (x,y,t).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} of the images {g.sup..sigma. (x,y,t).vertline.(x,y).epsilon..OMEGA.} corresponding to every parameter value .sigma..epsilon.[1, .varies.) and taken at a fixed time moment t.epsilon.T will be called the "parametric image" taken at the time moment t.

Each component g.sup..sigma. (x,y,t), g.epsilon.G of the parametric time-varying image function g.sup..sigma. (x,y,t) is infinitely differentiable everywhere in the domain .OMEGA.xT, and its partial derivatives with respect to the variables x, y, t can be obtained as observations of the combined generalized initial image function .GAMMA..sub.g on the partial derivatives with respect to the parameters x, y, t of the testing function .PSI..sub.x,y,t.sup..sigma. .epsilon..PHI.(R.sup.3) specified by the relation (2-4). For example, the components g.sub.x.sup..sigma. (x,y,t), g.sub.y.sup..sigma. (x,y,t), g.epsilon.G of the first-order partial derivatives g.sub.x.sup..sigma. (x,y,t).ident.{g.sub.x.sup..sigma. (x,y,t).vertline.g.epsilon.G}, g.sub.y.sup..sigma. (x,y,t).ident.{g.sub.y.sup..sigma. (x,y,t).vertline.g.epsilon.G} of the parametric time-varying image function g.sup..sigma. (x,y,t) are given as the observations ##EQU9## of the combined generalized initial image function .GAMMA..sub.g on the testing functions ##EQU10## for every (x,y,t).epsilon..OMEGA.xT, .sigma..epsilon.[1, .varies.), where ##EQU11## and .PSI..sub.x (x,y,t), .PSI..sub.y (x,y,t) are the partial derivatives of the function .PSI.(x,y,t) with respect to the variables x, y.

Let t-.DELTA.t.sup.-, t, t+.DELTA.t.sup.+ .epsilon.R be three consecutive instants of time (here .DELTA.t.sup.- and .DELTA.t.sup.+ are positive time increments), and let (x(t-.DELTA.t.sup.-), y(t-.DELTA.t.sup.-)).epsilon.R.sup.2 be a projection of some point in the scene taken at the time t-.DELTA.t.sup.-, while (x(t+.DELTA.t.sup.+),y(t+.DELTA.t.sup.+)).epsilon.R.sup.2 is a projection of the same point taken at the time t+.DELTA.t.sup.+. The vector (.mu..sub.A (x,y,t,.DELTA.t.sup.-,.DELTA.t.sup.+),.nu..sub.A (x,y,t,.DELTA.t.sup.-,.DELTA.t.sup.+)) defined by the relations:

.mu..sub.A (x,y,t,.DELTA.t.sup.-,.DELTA.t.sup.+)=(x(t+.DELTA.t.sup.+)-x(t-t.sup.-))/(.DELTA.t.sup.- +.DELTA.t.sup.+), (2-7) .nu..sub.A (x,y,t,.DELTA.t.sup.-,.DELTA.t.sup.+)=(y(t+.DELTA.t.sup.-)-y(t-.DELTA.t.sup.-))/(.DELTA.t.sup.- +.DELTA.t.sup.+), (2-8) where x=(.DELTA.t.sup.- x(t+.DELTA.t.sup.+)+.DELTA.t.sup.+ x(t-.DELTA.t.sup.-))/(.DELTA.t.sup.- +.DELTA.t.sup.+), (2-9) y=(.DELTA.t.sup.- y(t+.DELTA.t.sup.+)+.DELTA.t.sup.+ y(t-.DELTA.t.sup.-))/(.DELTA.t.sup.- +.DELTA.t.sup.+), (2-10) will be called the "average velocity vector" of the scene projection at the point (x,y).epsilon.R.sup.2 and the time t.epsilon.R based on the scene projections at the time t+.DELTA.t.sup.+ and time t-.DELTA.t.sup.-. The vector (.mu..sub.A (x,y,t,.DELTA.t.sup.-,.DELTA.t.sup.+),.nu..sub.A (x,y,t,.DELTA.t.sup.-,.DELTA.t.sup.+)) can be defined only at those points (x,y) from the set R.sup.2 that are projections of the points in the scene visible at the time t-.DELTA.t.sup.- as well as the time t+.DELTA.t.sup.+. The collection of such image points will be denoted by the symbol .PSI.(t,.DELTA.t.sup.-,.DELTA.t.sup.+).

The "instantaneous velocity vector" (.mu..sub.l (x,y,t),.nu..sub.l (x,y,t)) of the scene projection at the point (x,y).epsilon.R.sup.2 and time t.epsilon.R will now be defined. Let .DELTA.t.sup.-, .DELTA.t.sup.+ be positive time increments, and let W(x,y,t,.DELTA.t.sup.-,.DELTA.t.sup.+) be the set which is equal to a single-element set containing the average velocity vector (.mu..sub.A (x,y,t,.DELTA.t.sup.-,.DELTA.t.sup.+), .nu..sub.A (x,y,t,.DELTA.t.sup.-,.DELTA.t.sup.+)) if (x,y).epsilon..OMEGA.(t,.DELTA.t.sup.-,.DELTA.t.sup.+), and is equal to an empty set .O slashed. otherwise. Then for every positive time increment .DELTA.t<.delta.t (here .delta.t is a positive real constant), the set W(x,y,t.DELTA.t) will be defined from the relation ##EQU12## while the set W(x,y,t) will be defined as ##EQU13## where the symbol W(x,y,t,.DELTA.t) means the topological closure of the set W(x,y,t,.DELTA.t). It will be assumed that W(x,y,t) are single-element sets for almost every (x,y,t).epsilon.R.sup.3. This assumption means that the subset of the set R.sup.3 containing the points (x,y,t), where the sets W(x,y,t) are not single-element sets, has a Lebesgue measure equal to zero. For every point (x,y,t).epsilon.R.sup.3 such that the set W(x,y,t) consists of a single element, that element is selected as the value of the instantaneous velocity vector (.mu..sub.l (x,y,t),.nu..sub.l (x,y,t)); otherwise a zero is chosen for the value of the instantaneous velocity vector (.mu..sub.l (x,y,t),.nu..sub.l (x,y,t)).

Assume that the functions .mu..sub.l (x,y,t),.nu..sub.l (x,y,t), (x,y,t).epsilon.R.sup.3 are locally integrable with respect to the Lebesgue measure dx dy dt in R.sup.3 and thereby can be used to form continuous linear functional U.sub.l, V.sub.l (generalized functions) defined on the locally convex linear topological space .PHI.(R.sup.3). The values of the generalized functions U.sub.l, V.sub.l at the testing function .phi..epsilon..PHI.(R.sup.3) are given by the relations ##EQU14##

Let (x,y) be an image point from .OMEGA., t a moment of time from T, and .sigma. a parameter value from [1, .varies.); then the "velocity vector" (.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)) of the image corresponding to the parameter value .sigma. at the point (x,y) and time t is defined from the relations ##EQU15## where .PSI. is the measurement function from the set .PHI.(R.sup.3). The collection of the velocity vectors {(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.} corresponding to the moment of time t.epsilon.T and to the parameter value .sigma..epsilon.[1, .varies.) will be called the "velocity vector field" of the image {g.sup..sigma. (x,y,t).vertline.(x,y).epsilon..OMEGA.} corresponding to the parameter value .sigma. and taken at the time moment t. The collection {(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} of the velocity vector fields {(.mu..sup..sigma. (x,y,t),.mu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.} of the images {g.sup..sigma. (x,y,t).vertline.(x,y).epsilon..OMEGA.} corresponding to every parameter value .sigma..epsilon.[1, .varies.) and taken at a fixed time moment t.epsilon.T will be called the "parametric velocity vector field" of the parametric image {g.sup..sigma. (x,y,t).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} taken at the time moment t.

Let t.sub.0 "<t.sub.1 "<. . . <t.sub.K" " be a finite increasing sequence of time moments from some time interval T, and let t be a time moment from the same time interval T. Then the estimation problem can be formulated as follows. Given a sequence of parametric images {g.sup..sigma. (x,y,t.sub.K ").vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)}K"/k=0 taken at the time moments t.sub.k ", k=0, 1, . . . , K", which will be called the "parametric time-varying image sequence", and given the time moment t, find the parametric velocity vector field {(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} which is an estimate of the parametric velocity vector field {(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} of the parametric image {g.sup..sigma. (x,y,t).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} taken at the time moment t.

Although changes in the scene giving rise to the parametric velocity vector field {(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} are reflected in the parametric time-varying image sequence {g.sup..sigma. (x,y,t.sub.K ").vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)}K"/k=0 the relation between them is not necessarily unique. The same parametric time-varying image sequence can be associated with different parametric velocity vector fields and, vice versa, the same parametric velocity vector field can result from different parametric time-varying image sequences.

The method of the present invention, then, includes determining an estimate of the parametric velocity vector field corresponding to a given parametric time-varying image sequence. In order for this determination to be possible, specific assumptions have to be made about the scene being imaged and about the imaging process itself. The assumptions we make are described in the following. Based on these assumptions constraints are imposed on the estimate of the parametric velocity vector field. The determination of the estimate is then reduced to solving the system of equations arising from such constraints for a given parametric time-varying image sequence.

Let {g.sup..sigma. (x,y,t.sub.K ").vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)}K"/k=0 be given parametric time-varying image sequence taken at the time moments t.sub.k ", k=0, 1, . . . , K", which form an increasing sequence within the time interval T, and let t be a given time moment within the time interval T. We consider the problem of finding an estimate {(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} of the parametric velocity vector field {(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} as a function of the parametric time-varying image sequence {g.sup..sigma. (x,y,t.sub.K ").vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)}K"/k=0.

In what follows, the estimate {(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} of the parametric velocity vector field {(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} is determined as the parametric vector field satisfying a set of constraints. These constraints are based on the following assumptions:

1. The scene to be imaged is assumed to have near-constant and near-uniform illumination, which means that the changes in the incident illumination of each surface patch in the scene are small across the space and across the time, and are mainly due to the orientation of such surface patch relative to the light source.

2. It is assumed that changes in the incident illumination of each surface patch in the scene across the time due to the changes in the orientation of that surface patch are sufficiently small to be neglected.

3. The radiance of each surface point in the scene is assumed to be nearly Lambertian, which means that locally its dependence on the viewing position is small and can be neglected.

4. As has been mentioned in the previous section, the irradiance .xi.(x,y,t), .xi..epsilon..XI. of each point (x,y) in the image plane R.sup.2 at a time t.epsilon.R is assumed to be roughly proportional to the radiance of the point in the scene projecting into the point (x,y) in the image plane at the time t, and the proportionality coefficient, to be independent of the location (x,y) and the time t.

5. The scene is assumed to be made out of opague objects. The criteria for selecting feature points (x.sub..mu.,y.sub..mu.,t.sub..mu.).epsilon.M and the criteria for assigning feature values .eta.(x.sub..mu.,y.sub..mu.,t.sub..mu.), .eta..epsilon.H to each feature point (x.sub..mu.,y.sub..mu.,t.sub..mu.).epsilon.M are assumed to represent the true properties of the objects and to be independent of the objects' spatial attitudes.

6. The velocities of neighboring points of the objects in the scene are assumed to be similar, on one hand, and to change slowly with respect to time, on the other hand. In other words, it is assumed that the parametric time-varying velocity vector field of the surface points of each object in the scene varies smoothly across the space and across the time.

Let II.sub.t be a subset of the space R.sup.2 defined by the relation

II.sub.t .ident.{(.DELTA.t.sup.-,.DELTA.t.sup.+).vertline.(.DELTA.t.sup.-,.DELTA.t.sup.+).epsilon.R.sup.2, .DELTA.t.sup.- +.DELTA.t.sup.+ >0, (t-.DELTA.t.sup.-), (t+.DELTA.t.sup.+).epsilon.{t.sub.k ")K"/k=0}, (3-1) let P.sub.t be a subset of the set II.sub.t, and let G.sub.t .ident.GxP.sub.t be the Cartesian product of the set G on the set P.sub.t. The elements from the set G.sub.t are denoted by the symbol g.sub.t, and the set G.sub.t is assumed to be a measurable space with some measure dg.sub.t defined on the family B(G.sub.t) of Borel subsets of the set G.sub.t. In practice, the set G.sub.t is finite, and the measure dg.sub.t is a point measure on the finite set G.sub.t, but a more general approach is used for the sake of uniformity of presentation.

Let g.sub.t .ident.(g,.DELTA.t.sup.-,.DELTA.t.sup.+).epsilon.G.sub.t, .sigma..epsilon.[1, .varies.), and let g.sub.t.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.) be a function defined by the relation

g.sub.t.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.).ident.(g.sup..sigma. (x+.DELTA.t.sup.+ .mu..sup..sigma. (x,y,t),y+.DELTA.t.sup.+ .nu..sup..sigma. (x,y,t),t+.DELTA.t.sup.+)-g.sup..sigma. (x-.DELTA.t.sup.- .mu..sup..sigma. (x,y,t),y-.DELTA.t.sup.- .nu..sup..sigma. (x,y,t),t-.DELTA.t.sup.-))/(.DELTA.t.sup.- +.DELTA.t.sup.+), (3-2) for every (x,y).epsilon..OMEGA., where (.mu..sup..sigma.,.nu..sup..sigma.).ident.(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)) is the velocity vector of the image corresponding to the parameter value .sigma. at the image point (x,y) and the time moment t. Here, to simplify the notation, the symbol g.sub.t is used as the part of the notation g.sub.t.sup..sigma. for the above defined function (3-2) in addition to its role as the index for such a function.

Let (x,y).epsilon..OMEGA. be an image point that is taken at the time t projection of some point in the scene that does not belong to the occluding boundaries. The occluding boundary of an object is defined as the points in the scene belonging to the portion of the object that projects to its silhouette. The assumptions 1-5 made at the beginning of this section imply that the absolute value of the function g.sub.t.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.) is small. Therefore it is natural to use the function

(g.sub.t.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.)).sup.2 (3-3) as a part of the functional whose minimum specifies the estimate (.mu..sup..sigma.,.nu..sup..sigma.).ident.(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)) of the velocity vector (.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)) of the image corresponding to the parameter value .sigma.. It is customary for the function (3-3) to be called the "optical flow constraint" corresponding to the parameter value .sigma. and the index g.sub.t [2, 12-14, 18-21].

Note that if the time increments .DELTA.t.sup.-, .DELTA.t.sup.+ approach zero, then for an integer n greater than or equal to 2 the nth-order partial derivative of the function g.sub.t.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.) with respect to the components .mu..sup..sigma. .ident..mu..sup..sigma. (x,y,t), .nu..sup..sigma. .ident..nu..sup..sigma. (x,y,t) of the estimate of the velocity vector approaches zero at a rate proportional to (.DELTA.t.sup.-).sup.n-1 +(.DELTA.t.sup.+).sup.n-1. This, in particular, implies that when the time increments .DELTA.t.sup.-, .DELTA.t.sup.+ approach zero the function g.sub.t.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.) approaches the function which is linear with respect to the components .mu..sup..sigma., .nu..sup..sigma. of the estimate of the velocity vector.

It can be easily verified that the estimate (.mu..sup..sigma. (x,y,t),.mu..sup..sigma. (x,y,t)) of the velocity vector on which the function (3-3) achieves its minimum satisfies the system of equations ##EQU16## where the function g.sub.t.mu..sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.) and g.sub.t.nu..sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.) are the first-order partial derivatives of the function g.sub.t.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.) with respect to the components .mu..sup..sigma. and .nu..sup..sigma. of the estimate of the velocity vector.

Let the vector (.DELTA..mu..sup..sigma. (x,y,t),.DELTA..nu..sup..sigma. (x,y,t)) be defined as in the relations

.mu..sup..sigma. (x,y,t)=.mu..sup..sigma. (x,y,t)+.DELTA..mu..sup..sigma.(x,y,t), (3-5) .nu..sup..sigma. (x,y,t)=.nu..sup..sigma. (x,y,t)+.DELTA..nu..sup..sigma. (x,y,t). (3-6) By expanding the left hand sides of the relation (3-4) into their Taylor series about the velocity vector (.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)) and by setting the second-order and higher partial derivatives of the function g.sub.t.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.) with respect to the components .mu..sup..sigma. .ident..mu..sup..sigma. (x,y,t) and .nu..sup..sigma. .ident..nu..sup..sigma. (x,y,t) of the velocity vector to zero, we obtain the system of linear equations, relative to the vector (.DELTA..mu..sup..sigma. (x,y,t),.DELTA..nu..sup..sigma. (x,y,t)), ##EQU17## where the functions g.sub.tu.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.) and g.sub.tv.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.) are the first-order partial derivatives of the function g.sub.t.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.) with respect to the components .mu..sup..sigma. and .nu..sup..sigma. of the velocity vector.

Let us first assume that the norm of the gradient

(g.sub.tu.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.),g.sub.tv.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.)) (3-8) of the function g.sub.t.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.), with respect to the components .mu..sup..sigma. and .nu..sup..sigma. of the velocity vector, is relatively large; then the system of equations (3-7) is equivalent to the relation g.sub.tu.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.).DELTA..mu..sup..sigma. (x,y,t)+g.sub.tv.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.).DELTA..nu..sup..sigma. (x,y,t)+g.sub.t.sup..sigma. (x,y,t,.mu..sup..sigma.,.nu..sup..sigma.)=0. (3-9) The last relation implies that the requirement for the function (3-3) to achieve a minimum on the estimate (.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)) of the velocity vector (.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)) specifies the component of the estimate in the direction of the gradient (3-8) and leaves the component orthogonal to the gradient (3-8) undetermined.

If the norm of the gradient (3-8) is relatively small; then the system of equations (3-7) becomes too weak, and the requirement for the function (3-3) to achieve a minimum on the estimate (.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)) of the velocity vector (.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)) imposes no constraint on the estimate.

The above discussion suggests that we need to impose additional constraints on the estimate {(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} of the parametric velocity vector field {(.mu..sup..sigma. (x,y,t),.nu..sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1, .varies.)} in order for the computation to be well defined.

Let S.ident.{(s.sub.x,s.sub.y).vertline.(s.sub.x,s.sub.y).epsilon.R.sup.2, x.sub.x.sup.2 +s.sub.y.sup.2 =1} be a unit circle in the two-dimensional Euclidean space R.sup.2, with each vector s.ident.(s.sub.x,s.sub.y).epsilon.S specifying the direction in the image plane. The set S is assumed to be a measurable space with some uniform measure ds defined on the family B(S) of Borel subsets of the set S.

Let s.tbd.(s.sub.x,s.sub.y).epsilon.S, and let (s,.gradient.u.sup..sigma. (x,y,t)), (s,.gradient.v.sup..sigma. (x,y,t)) be the functions defined by the relations ##EQU18## for every (x,y).epsilon..OMEGA.. Here the symbol .gradient. stands for the gradient of a function with respect to the variables x,y.

Let s.tbd.(s.sub.x,s.sub.y).epsilon.S, and let (x,y).epsilon..OMEGA. be an image point that is taken at a time t projection of some point in the scene that does not cross the occluding boundary in the direction s, meaning that for all sufficiently small positive constants .omega., the following two points in the scene both belong to the same object and are sufficiently close to each other: the first is the one projecting to the image point (s+.omega.s.sub.x,y+.omega.s.sub.y) at the time t, whereas the second is the one projecting to the image point (x,y) at the time t. The assumption 6 made at the beginning of this section implies that the absolute values of the functions (s,.gradient.u.sup..sigma. (x,y,t)) and (s,.gradient.v.sup..sigma. (x,y,t)) are small. Therefore it is natural to use the functions

(s,.gradient.u.sup..sigma. (x,y,t)).sup.2, (3-12) (s,.gradient.v.sup..sigma. (x,y,t)).sup.2 (3-13) as a part of the functional whose minimum specifies the estimate (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) of the velocity vector (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) of the image corresponding to the parameter value .sigma.. It is customary for the functions (3-12), (3-13) to be called the "smoothness constraints" in the direction s [2, 12-14, 18-21].

Finally, to make the process of computing the estimate {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} well defined, we require that for every image point (x,y).epsilon..OMEGA. and for every parameter value .sigma..epsilon.[1,,.infin.) where the optical flow constraints and the directional smoothness constraints impose no restrictions the difference between the unknown estimate (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) of the velocity vector (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) and its initial estimate (u.sub.0.sup..sigma. (x,y,t),v.sub.0.sup..sigma. (x,y,t)) be small. This can be accomplished by including the functions

(u.sup..sigma. (x,y,t)-u.sub.0.sup..sigma. (x,y,t)).sup.2, (3-14) (v.sup..sigma. (x,y,t)-v.sub.0.sup..sigma. (x,y,t)).sup.2 (3-15) in the functional whose minimum specifies the estimate (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) of the velocity vector (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) of the image corresponding to the parameter value .sigma. in addition to the functions (3--3), (3-12), and (3-13) specifying the optical flow and the directional smoothness constraints, provided that the weights associated with the functions (3-14) and (3-15) are small relative to the weights associated with the functions (3--3), (3-12), and (3-13). The initial estimate (u.sub.0.sup..sigma. (x,y,t),v.sub.0.sup..sigma. (x,y,t)) appearing in the relations (3-14), (3-15) is defined towards the end of the next section. The constraints imposed on the estimate (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) of the velocity vector (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) at the image point (x,y).epsilon..OMEGA. and the time moment t by the above requirements will be called the "regularization constraints".

As has been stated earlier, the optical flow and the directional smoothness constraints are not necessarily valid at the points near the occluding boundaries, even when the assumptions described at the beginning of this section are observed. The method of computing the estimate of the parametric velocity vector field of the present invention resolves the above difficulties by adjusting the weight associated with each constraint in such a way that it becomes small whenever the constraint is not valid. In the following, it is described how the functions (3--3), (3-12), (3-13), (3-14), and (3-15), specifying the optical flow, directional smoothness, and regularization constraints, are combined into the functional of the estimate of the parametric velocity vector field. The estimate is then computed by solving the system of nonlinear equations arising from a certain optimality criterion related to such functional.

The following describes the system of nonlinear equations used in the method of the present invention in its estimation of velocity vector fields. In the method of the invention, the estimate {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} is determined as the parametric vector field on which a weighted average of the optical flow, directional smoothness, and regularization constraints is minimized. The weight functions are chosen in such a way that at every image point (x,y).epsilon..OMEGA. which is near the occluding boundary the following requirements are met: each weight function associated with an optical flow constraint becomes small whenever the optical flow constraint is not satisfied, and each weight function associated with a smoothness constraint corresponding to a direction in which an occluding boundary is crossed becomes small whenever the directional smoothness constraint is not satisfied. The criteria for the presence of the occluding boundary near a given image point, and the criteria for the occluding boundary being crossed near a given image point in a given direction, can be defined most effectively in terms of the values of the unknown estimate of the parametric velocity vector field. Therefore, the values of the unknown estimate of the parametric velocity vector field have to be used implicitly in each of the weight functions. On the other hand, each weight function has to be treated as if it were independent of the values of the unknown estimate of the parametric velocity vector field, because it only specifies a relative significance of the corresponding constraint as a part of the weighted average and not the constraint itself. To overcome these difficulties, two copies of the unknown estimate of the parametric velocity vector field are introduced: the invariable one, and the variable one. In the weight functions the values of the invariable copy of the unknown estimate of the parametric velocity vector field are used, whereas in the constraint functions the values of the variable copy of the unknown estimate of the parametric velocity vector field are used. Then a more general variational principle as opposed to an energy minimization principle is applied to the weighted average of the optical flow, directional smoothness, and regularization constraints to derive the system of nonlinear equations of the unknown estimate of the parametric velocity vector field. This variational principle can be described as follows.

Let (u,v).tbd.{(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} be an invariable copy of the estimate of the parametric velocity vector field, let (u,v).tbd.{(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} be a variable copy of the estimate of the parametric velocity vector field, and let the parametric vector field (.increment.u,.increment.v).tbd.{(.increment.u.sup..sigma. (x,y,t),.increment.v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} be defined by the relation

(.increment.u.sup..sigma. (x,y,t),.increment.v.sup..sigma. (x,y,t))=(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t))-(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)), (x,y).epsilon..OMEGA., .sigma..epsilon.[1,.infin.) (4-1) Assume that the derivative of the parametric vector field (.increment.u,.increment.v) in the direction normal to the boundary .differential..OMEGA. of the image .OMEGA. is equal to zero for every parameter value .sigma..epsilon.[1,.infin.) and for every image point (x,y) belonging to the image boundary .differential..OMEGA.. This will impose restrictions on the types of the variations that can be applied to the variable copy of the unknown estimate of the parametric velocity vector field.

Let f(u,v,u,v) be a functional of the parametric vector fields (u,v), (u,v) defined as a weighted average of functions (3--3), (3-12), (3-13), (3-14), and (3-15), specifying the optical flow, directional smoothness, and regularization constraints, respectively, by the following relation: ##EQU19## Here: (i) .alpha..sub.g.sbsb.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,.parallel..gradient.u.sup..sigma. .parallel.,.parallel..gradient.v.sup..sigma. .parallel.), (x,y).epsilon..OMEGA., .sigma..epsilon.[1,.infin.), g.sub.t .epsilon.G.sub.t is a weight associated with an optical flow constraint corresponding to the parameter value .sigma. and the index g.sub.t, which is a functional of the independent variables x, y, t and of the functions u.sup..sigma. .tbd.u.sup..sigma. (x,y,t), v.sup..sigma. .tbd.v.sup..sigma. (x,y,t), .parallel..gradient.u.sup..sigma. .parallel..tbd..parallel..gradient.u.sup..sigma. (x,y,t).parallel., .parallel..gradient.v.sup..sigma. .parallel..tbd..parallel..gradient.v.sup..sigma. (x,y,t).parallel. of these variables;

(ii) .beta..sub.s.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,(s,.gradient.u.sup..sigma.),(s,.gradient.v.sup..sigma.)), (x,y).epsilon..OMEGA., .sigma..epsilon.[1,.infin.), s.epsilon.S is a weight associated with a smoothness constraint in the direction s, which is a function of the independent variables x, y, t and of the functions u.sup..sigma. .tbd.u.sup..sigma. (x,y,t), v.sup..sigma. .tbd.v.sup..sigma. (x,y,t), (s,.gradient.u.sup..sigma.).tbd.(s,.gradient.u.sup..sigma. (x,y,t)), (s,.gradient.v.sup..sigma.).tbd.(s,.gradient.v.sup..sigma. (x,y,t)) of these variables; (iii) .gamma..sub.0 is a positive constant specifying a weight for the regularization constraint; (iv) the functions .parallel..gradient.u.sup..sigma. (x,y,t).parallel., .parallel..gradient.v.sup..sigma. (x,y,t).parallel., (x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.) are the norms of the vector valued functions .gradient.u.sup..sigma. (x,y,t), .gradient.v.sup..sigma. (x,y,t), which are the gradients with respect to the variables x,y of the functions u.sup..sigma. (x,y,t), v.sup..sigma. (x,y,t), respectively, given by the relations ##EQU20## (v) the functions (s,.gradient.u.sup..sigma. (x,y,t)), (s,.gradient.v.sup..sigma. (x,y,t)), (s,.gradient.u.sup..sigma. (x,y,t)), (s,.gradient.v.sup..sigma. (x,y,t)), (x,y).epsilon..OMEGA., .sigma..epsilon.[1,.infin.), s.epsilon.S are defined as in the relations (3-10), (3-11); (vi) the function g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.), (x,y).epsilon..OMEGA., .sigma..epsilon.[1,.infin.), g.sub.t .epsilon.G.sub.t is defined as in the relation (3-2); and (vii) the components u.sub.0.sup..sigma. (x,y,t), v.sub.0.sup..sigma. (x,y,t) of the initial estimate of the velocity vector corresponding to the parameter value .sigma..epsilon.[1,.infin.) at the image point (x,y).epsilon..OMEGA. are defined later.

The estimate of the parametric velocity vector field is then defined as the parametric vector field (u,v), on which the functional f(u,v,u,v), considered as the function of the parametric vector field (u,v) and depending on the parametric vector field (u,v) as on the parameters, achieves a local minimum when the value of the parametric vector field (u,v) is identically equal to the value of the parametric vector field (u,v). Taking into account the relation (4-1) the functional f(u,v,u,v) can be expressed in the form f(u,v,u+.increment.u,v+.increment.v). The parametric vector field (.increment.u,.increment.v) specifies a perturbation to the parametric vector field (u,v), and the functional f(u,v,u+.increment.u,v+.increment.v) assigns the cost to each choice of the parametric vector field (u,v) and its perturbation (.increment.u,.increment.v). Then the estimate of the parametric velocity vector field is the parametric vector field (u,v), for which a locally minimal cost is achieved when the perturbation (.increment.u,.increment.v) is identically equal to zero. Using the calculus of variations with respect to the parametric vector field (.increment.u,.increment.v) applied to the functional f(u,v,u+.increment.u,v+.increment.v), the above defined estimate (u,v) of the parametric velocity vector field is a solution of the system of equations ##EQU21## where the functions g.sub.tu.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) and g.sub.tv.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) are the first-order partial derivatives of the function g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) with respect to the components u.sup..sigma. .tbd.u.sup..sigma. (x,y,t) and v.sup..sigma. .tbd.v.sup..sigma. (x,y,t) of the estimate of the velocity vector.

To select the weight functions associated with the optical flow, directional smoothness, and regularization constraints, consider the functional f(u,v,u,v) defined by the relation (4-2) under the conditions that the parametric vector field (u,v) be identically equal to the parametric vector field (u,v), and that both be equal to the estimate {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1,.infin.)} of the parametric velocity vector field obtained as the solution of the system of nonlinear equations (4-5). In this case the functional f(u,v,u,v) takes the form ##EQU22##

Let {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1,.infin.)} be an estimate of the parametric velocity vector field, which is obtained as the solution of the system of nonlinear equations (4-5). Select each of the weight functions .alpha..sub.g.sbsb.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,.parallel..gradient.u.sup..sigma. .parallel.,.parallel..gradient.v.sup..sigma. .parallel.), g.sub.t .epsilon.G.sub.t in such a way that the contribution of the optical flow constraint g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) to the functional (4-6) becomes small for every image point (x,y).epsilon..OMEGA. located near the occluding boundary where the optical flow constraint g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) is not satisfied. Similarly, each of the weight functions .beta..sub.s.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,(s,.gradient.u.sup..sigma.), (s,.gradient.v.sup..sigma.)), s.epsilon.S will be chosen so that the contributions of the directional smoothness constraints (s,.gradient.u.sup..sigma. (x,y,t)) and (s,.gradient.v.sup..sigma. (x,y,t)) to the functional (4-6) become small for every image point (x,y).epsilon..OMEGA. located near the occluding boundary where the occluding boundary is crossed in the direction s and one of the directional smoothness constraints (s,.gradient.u.sup..sigma. (x,y,t)) and (s,.gradient.v.sup..sigma. (x,y,t)) is not satisfied. Let the image point (x,y).epsilon..OMEGA. be near the occluding boundary; then the following two events are likely to happen:

1. For some small positive time increments .increment.t.sup.- and .increment.t.sup.+ the point in the scene projecting into the image point (x,y) at the time t either is visible at the time (t-.increment.t.sup.-) and is invisible at the time (t+.increment.t.sup.+), or it is visible at the time (t+.increment.t.sup.+) and is invisible at the time (t-.increment.t.sup.-). Let (x.sub.1,y.sub.1) and x.sub.2,y.sub.2) be projections of such a point in the scene into the image plane taken at the time moments (t-.increment.t.sup.-) and (t+.increment.t.sup.+), respectively. If the radiance of the point in the scene projecting into the image point (x.sub.1,y.sub.1) at the time (t-.increment.t.sup.-) is significantly different from the radiance of the point in the scene projecting into the image point (x.sub.2,y.sub.2) at the time (t+.increment.t.sup.+), and if the radiance undergoes some changes in the neighborhoods of the above points in the scene, then the following cases are likely to be encountered:

1.1. The point in the scene projecting into the image point (x-.increment.t.sup.- u.sup..sigma. (x,y,t),y-.increment.t.sup.- v.sup..sigma. (x,y,t)) at the time (t-.increment.t.sup.-), and the point in the scene projecting into the image point (x+.increment.t.sup.+ u.sup..sigma. (x,y,t),y+.increment.t.sup.+ v.sup..sigma. (x,y,t)) at the time (t+.increment.t.sup.+), both belong to the occluded object, on one hand, and have different radiances, on the other hand. The later may happen if the above point occupy distinct locations on the object of the scene which, in turn, may be caused by the presence of the occluding boundary. Under these conditions the absolute values of some of the functions g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.), g.sub.t .epsilon.G.sub.t become large.

1.2. The point in the scene projecting into the image point (x-.increment.t.sup.- u.sup..sigma. (x,y,t),y-.increment.t.sup.- v.sup..sigma. (x,y,t)) at the time (t-.increment.t.sup.-) and the point in the scene projecting into the image point (x+.increment.t.sup.+ u.sup..sigma. (x,y,t),y+.increment.t.sup.+ v.sup..sigma. (x,y,t)) at the time (t+.increment.t.sup.+) both belong either to the occluding object of to the occluded one. In addition to this, for some not-too-large vector (.increment.u.sup..sigma. (x,y,t),.increment.v.sup..sigma. (x,y,t)), one of the following two points in the scene belongs to the occluding object, while the other belongs to the occluded one: the first point is the one projecting into the image point (x-.increment.t.sup.- (u.sup..sigma. (x,y,t)+.increment.u.sup..sigma. (x,y,t)),y-.increment.t.sup.- (v.sup..sigma. (x,y,t)+.increment.v.sup..sigma. (x,y,t))) at the time (t-.increment.t.sup.-), and the second point is the one projecting into the image point (x+.increment.t.sup.+ (u.sup..sigma. (x,y,t)+.increment.u.sup..sigma. (x,y,t)),y+.increment.t.sup.+ (v.sup..sigma. (x,y,t)+.increment.v.sup..sigma. (x,y,t))) at the time (t+.increment.t.sup.+). Likewise, for the same vector (.increment.u.sup..sigma. (x,y,t),.increment.v.sup..sigma. (x,y,t)), one of the following two points in the scene belongs to the occluding object, while the other belongs to the occluded one: the first point is the one projecting into the image point (x-.increment.t.sup.- (u.sup..sigma. (x,y,t)-.increment.u.sup..sigma. (x,y,t)),y-.increment.t.sup.- (v.sup..sigma. (x,y,t)-.increment.v.sup..sigma. (x,y,t))) at the time (t-.increment.t.sup.-), and the second point is the one projecting into the image point (x+.increment.t.sup.+ (u.sup..sigma. (x,y,t)-.increment.u.sup..sigma. (x,y,t)),y+.increment.t.sup.+ (v.sup..sigma. (x,y,t)-.increment.v.sup..sigma. (x,y,t))) at the time (t+.increment.t.sup.+) (see FIG. 1). In such a case the values of the functions (g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2, g.sub.t .epsilon.G.sub.t become relatively small, while the values of some of the functions (g.sub.t.sup..sigma. (x,y,t,u.sup..sigma. +.increment.u.sup..sigma.,v.sup..sigma. +.increment.v.sup..sigma.)).sup.2, g.sub.t .epsilon.G.sub.t, as well as the values of some of the functions (g.sub.t.sup..sigma. (x,y,t,u.sup..sigma. -.increment.u.sup..sigma.,v.sup..sigma. -.increment.v.sup..sigma.)).sup.2, g.sub.t .epsilon.G.sub.t, become relatively large.

This implies that the values of some of the functions

(g.sub.t.sup..sigma. (x,y,t,u.sup..sigma. +.increment.u.sup..sigma.,v.sup..sigma. +.increment.v.sup..sigma.)).sup.2 +(g.sub.t.sup..sigma. (x,y,t,u.sup..sigma. -.increment.u.sup..sigma.,v.sup..sigma. -.increment.v.sup..sigma.)).sup.2 -2(g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2, (4-7) g.sub.t .epsilon.G.sub.t becomes large. If the second-order and higher partial derivatives of the function g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) with respect to the components u.sup..sigma. .tbd.u.sup..sigma. (x,y,t) and v.sup..sigma. .tbd.v.sup..sigma. (x,y,t) of the estimate of the velocity vector are ignored, each of the functions (4-7) can be expressed in the form 2(.increment.u.sup..sigma. (x,y,t)g.sub.tu.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)+.increment.v.sup..sigma. (x,y,t)g.sub.tv.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2 (4-8) where the functions g.sub.tu.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) and g.sub.tv.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) are the first-order partial derivatives of the function g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) with respect to the components u.sup..sigma. and v.sup..sigma. of the estimate of the velocity vector. For the vector s.tbd.(s.sub.x,s.sub.y), s.epsilon.S that is collinear to the vector (.increment.u.sup..sigma. (x,y,t),.increment.v.sup..sigma. (x,y,t)), the above observations imply that the absolute values of some of the functions (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)), g.sub.t .epsilon.G.sub.t, defined by the relation (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.))=s.sub.x g.sub.tu.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)+s.sub.y g.sub.tv.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) (4-9) become large.

2. The estimate {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma..epsilon.[1,.infin.)} of the parametric velocity vector field changes rapidly in the vicinity of the image point (x,y) (see FIG. 2), which implies the following:

2.1. The values of the functions .parallel..gradient.u.sup..sigma. (x,y,t).parallel. and .parallel..gradient.v.sup..sigma. (x,y,t).parallel. become large.

2.2. The absolute values of the functions (s,.gradient.u.sup..sigma. (x,y,t)) and (s,.gradient.v.sup..sigma. (x,y,t)), s.epsilon.S become large, provided that the image point (x,y) crosses the occluding boundary in the direction s.

If the conditions on the scene being imaged and on the imaging process itself described at the beginning of the previous section are satisfied, then in the event of the absolute value of the function g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) being large for some g.sub.t .epsilon.G.sub.t, or of the values of the functions .parallel..gradient.u.sup..sigma. (x,y,t).parallel., .parallel..gradient.v.sup..sigma. (x,y,t).parallel. being large in addition to the absolute value of the function g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) being large for some g.sub.t .epsilon.G.sub.t, the image point (x,y) is likely to be near the occluding boundary, while in the event of the absolute values of the functions (s,.gradient.u.sup..sigma. (x,y,t)), (s,.gradient.v.sup..sigma. (x,y,t)), s.epsilon.S being large, or of the absolute values of some of the functions (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)), g.sub.t .epsilon.G.sub.t being large in addition to the absolute values of the functions (s,.gradient.u.sup..sigma. (x,y,t)), (s,.gradient.v.sup..sigma. (x,y,t)), s.epsilon.S being large, the image point (x,y) is likely to cross the occluding boundary in the direction s. Note that in the case of the functions .parallel..gradient.u.sup..sigma. (x,y,t).parallel., .parallel..gradient.v.sup..sigma. (x,y,t).parallel. being large, the image point (x,y) does not necessarily lie near the occluding boundary. It may, for example, lie on the object of the scene whose projection onto the image plane undergoes a rotation or a local deformation. Also note that in the case of some of the functions (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)), g.sub.t .epsilon.G.sub.t being large, the image point (x,y) does not necessarily cross the occluding boundary in the direction s. It may, for example, cross the radiance boundary arising from a texture or a sudden change in the illumination.

These observations suggest the following: each of the weight functions .alpha..sub.g.sbsb.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,.parallel..gradient.u.sup..sigma. .parallel.,.parallel..gradient.v.sup..sigma. .parallel.), g.sub.t .epsilon.G.sub.t should be a steadily decreasing function relative to the absolute value of the function g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) and to the values of the functions .parallel..gradient.u.sup..sigma. (x,y,t).parallel., .parallel..gradient.v.sup..sigma. (x,y,t).parallel. multiplied by the absolute value of the function g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.); and each of the weight functions .beta..sub.s.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,(s,.gradient.u.sup..sigma.),(s,.gradient.v.sup..sigma.)), s.epsilon.S should be a steadily decreasing function relative to the absolute values of the functions (s,.gradient.u.sup..sigma. (x,y,t)), (s,.gradient.v.sup..sigma. (x,y,t)) and to the absolute values of some of the functions (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)), g.sub.t .epsilon.G.sub.t multiplied by some combination of the absolute values of the functions (s,.gradient.u.sup..sigma. (x,y,t)), (s,.gradient.v.sup..sigma. (x,y,t)).

It is natural to require each of the weight functions corresponding to either an optical flow constraint or a directional smoothness constraint to be selected in such a way that the product of the weight function and the corresponding constraint function appearing as the summand in the functional (4-6) becomes larger whenever the constraint function becomes larger. This makes the solution of the system of nonlinear equations (4-5) determining the estimate of the parametric velocity vector field more robust with respect to the choice of the initial estimate of the parametric velocity vector field and with respect to the noise. To be more specific, these requirements consist of the following.

1. For every g.sub.t .epsilon.G.sub.t the function defined by the relation:

.alpha..sub.g.sbsb.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,.parallel..gradient.u.sup..sigma. .parallel.,.parallel..gradient.v.sup..sigma. .parallel.)(g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2,(4-10) and appearing as the summand in the functional (4-6) is required to be steadily increasing relative to the function (g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2. In other words, the function (4-10) is required to increase whenever the function (g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2 is increased.

2. For every s.epsilon.S the function defined by the relation:

.beta..sub.s.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,(s,.gradient.u.sup..sigma.),(s,.gradient.v.sup..sigma.))((s,.gradient.u.sup..sigma. (x,y,t)).sup.2 +(s,.gradient.v.sup..sigma. (x,y,t)).sup.2) (4-11) and appearing as the summand in the functional (4-6) is required to be steadily increasing relative to the function ((s,.gradient.u.sup..sigma. (x,y,t)).sup.2 +(s,.gradient.v.sup..sigma. (x,y,t)).sup.2). In other words, the function (4-11) is required to increase whenever the function ((s,.gradient.u.sup..sigma. (x,y,t)).sup.2 +(s,.gradient.v.sup..sigma. (x,y,t)).sup.2) is increased.

The weight functions

.alpha..sub.g.sbsb.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,.parallel..gradient.u.sup..sigma. .parallel.,.parallel..gradient.v.sup..sigma. .parallel.)=(r.sup.2 +(p.sup.2 +q.sup.2 (.parallel..gradient.u.sup..sigma. (x,y,t).parallel..sup.2 +.parallel..gradient.v.sup..sigma. (x,y,t).parallel..sup.2))(g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2).sup.-1, (4-12) for every g.sub.t .epsilon.G.sub.t, where r.sup.2, p.sup.2, q.sup.2 are non-negative real constants, and .beta..sub.s.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,(s,.gradient.u.sup..sigma.),(s,.gradient.v.sup..sigma.))=(a.sup.2 +(c.sup.2 +b.sup.2 (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.).sup.2)((s,.gradient.u.sup..sigma. (x,y,t)).sup.2 +(s,.gradient.v.sup..sigma. (x,y,t)).sup.2)).sup.-1, (4-13) for every s.epsilon.S, where the expression b.sup.2 (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2 is defined as in the relation ##EQU23## and a.sup.2, c.sup.2, b.sub.g.sbsb.t.sup.2, g.sub.t .epsilon.G.sub.t are non-negative real constants, comply with the above requirements and with the requirements stated earlier in this section.

The approach of obtaining the estimate of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} by solving the system of nonlinear equations (4-5) for each parameter value .sigma. independently may encounter difficulty relating to the fact that the system of equations (4-5) may, in general, have many solutions. To overcome this difficulty, the values of the estimates of the velocity vector field corresponding to the different values of the parameter .sigma. are tied together by imposing the following restriction: the estimate of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} is required to be a continuously differentiable function with respect to the parameter .sigma. for every image point (x,y).epsilon..OMEGA. and for every parameter value .sigma..epsilon.[1,.infin.). The additional restrictions on the estimate of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} are imposed in the form of the boundary conditions as follows:

1. The estimate of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} is required to converge to the vector field that is identically equal to a given vector constant (u.sup..infin.,v.sup..infin.) when the parameter .sigma. converges to .infin., and such convergence is required to be uniform over the set of the image points .OMEGA..

2. The derivative of the estimate of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} in the direction normal to the boundary .differential..OMEGA. of the image .OMEGA. is required to be equal to zero for every image point (x,y) belonging to the image boundary .differential..OMEGA. and for every parameter value .sigma..epsilon.[1,.infin.).

Let .sigma..epsilon.[1,.infin.) be a given parameter value, and let .increment..sigma. be an infinitesimal positive real value. Taking into consideration the above restrictions, the initial estimate of the velocity vector field {(u.sub.0.sup..sigma. (x,y,t),v.sub.0.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.} corresponding to the parameter value .sigma. can be defined as follows:

(u.sub.0.sup..sigma. (x,y,t),v.sub.0.sup..sigma. (x,y,t)).tbd.(u.sup..sigma.+.increment..sigma. (x,y,t),v.sup..sigma.+.increment..sigma. (x,y,t)), (x,y).epsilon..OMEGA..(4-15) Then the regularization constraints (3-14) and (3-15) can, respectively, be approximated by the relations ##EQU24## for every image point (x,y).epsilon..OMEGA., where .gamma. is a positive real constant given by the relation

By substituting the regularization constraints .gamma..sub.0 (u.sup..sigma. (x,y,t)-u.sub.0.sup..sigma. (x,y,t)) and .gamma..sub.0 (v.sup..sigma. (x,y,t)-v.sub.0.sup..sigma. (x,y,t)) appearing in the system of nonlinear equations (4-5) with the constraints ##EQU25## respectively, and by defining the weight functions .alpha..sub.g.sbsb.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,.parallel..gradient.u.sup..sigma. .parallel.,.parallel..gradient.v.sup..sigma. .parallel.), g.sub.t .epsilon.G.sub.t and .beta..sub.s.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,(s,.gradient.u.sup..sigma.),(s,.gradient.v.sup..sigma.)), s.epsilon.S according to the relations (4-12) and (4-13), respectively, the system (4-5) becomes the following system of equations relative to the unknown estimate of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)}: ##EQU26##

Under conditions that the regularization constraints 0.5.gamma..sub.0 (u.sup..sigma. (x,y,t)-u.sub.0.sup..sigma. (x,y,t)).sup.2 and 0.5.gamma..sub.0 (v.sup..sigma. (x,y,t)-v.sub.0.sup..sigma. (x,y,t)).sup.2 appearing in the functional f(u,v,u,v) defined by the relation (4-2) be replaced with the constraints ##EQU27## respectively, and provided that the weight functions .alpha..sub.g.sbsb.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,.parallel..gradient.u.sup..sigma. .parallel.,.parallel..gradient.v.sup..sigma. .parallel.), g.sub.t .epsilon.G.sub.t are defined according to the relation (4-12), whereas the weight functions .beta..sub.s.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,(s,.gradient.u.sup..sigma.),(s,.gradient.v.sup..sigma.)), s.epsilon.S are defined according to the relation (4-13), the functional f(u,v,u,v) of the parametric vector field (u,v).tbd.{(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y.epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} and the parametric vector field (u,v).tbd.{(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} defined by the relation (4-2) can be replaced with the functional h(u,v,u,v) of the same parametric vector fields (u,v), (u,v) defined as in the relation ##EQU28## where the symbol ln stands for the natural logarithm function.

The above substitutions imply that the parametric vector field (u,v) satisfies the system of nonlinear equations (4-19) if and only if every variation of the functional h(u,v, u,v) with respect to the parametric vector field (u,v) is equal to zero whenever the parametric vector field (u,v) is identically equal to the parametric vector field (u,v).

Note that when the constant q.sup.2 approaches zero and when the increments between successive time moments of the time sequence {t.sub.k "}.sub.k= 0.sup.K"+0 uniformly approach zero the functional h(u,v,u,v) approaches the functional, which is independent of the parametric vector field (u,v) and the problem of finding the estimate of the parametric velocity vector field becomes the energy minimization problem.

Let (u,v).tbd.{(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} be the estimate of the parametric velocity vector field obtained as the solution of the system of nonlinear equations (4-19). We now examine the behavior of the system of nonlinear equations (4-19) near such a solution. First, the changes in the behavior of the optical flow and the directional smoothness constraints when the image point (x,y) goes from the position beyond the occluding boundary to the position near the occluding boundary are considered. Then, the dependence of the optical flow and the directional smoothness constraints on the parameters r.sup.2,p.sup.2,q.sup.2 specifying the weight functions .alpha..sub.g.sbsb.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,.parallel..gradient.u.sup..sigma. .parallel.,.parallel..gradient.v.sup..sigma. .parallel.), g.sub.t .epsilon.G.sub.t as in the relation (4-12) and on the parameters a.sup.2, c.sup.2,b.sub.g.sbsb.t.sup.2,g.sub.t .epsilon.G.sub.t specifying the weight functions .beta..sub.s.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,(s,.gradient.u.sup..sigma.),(s,.gradient.v.sup..sigma.)), s.epsilon.S as in the relation (4-13) are investigated.

Let (.increment.u,.increment.v).tbd.{(.increment.u.sup..sigma. (x,y,t),.increment.v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} be a parametric vector field defined by the relation (4-1), then in the neighborhood of the parametric vector field (u,v) the functional h(u,v,u+.increment.u,v+.increment.v) can be expanded into the quadratic form of the parametric vector field (.increment.u,.increment.v) as follows: ##EQU29##

The optical flow constraint g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) corresponding to the index g.sub.t .epsilon.G.sub.t is represented in the system of nonlinear equations (4-19) by the expression ##EQU30## With respect to this expression the following cases can be considered:

1. The image point (x,y) is beyond the occluding boundary. Then the value of the function (g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2 may become small. In this case the value of the function (4-22) would approach the value of the function ##EQU31## This, in turn, would make each of the two summands of the system of nonlinear equations (4-19) associated with the optical flow constraint g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) corresponding to the index g.sub.t behave as a linear function relative to the function g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) with the proportionality coefficient equal to the function g.sub.tu.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)/r.sup.2 in the case of the first summand and to the function g.sub.tv.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)/r.sup.2 in the case of the second summand.

2. The image point (x,y) approaches the occluding boundary. Then the value of the function (g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2 may increase. In the case when the value of the function (g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2 is less than the value of the function r.sup.2 /(p.sup.2 +q.sup.2 (.parallel..gradient.u.sup..sigma. (x,y,t).parallel..sup.2 +.parallel..gradient.v.sup..sigma. (x,y,t).parallel..sup.2)), the fact that the summand of the quadratic form (4-21) associated with the optical flow constraint g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) is positive implies that the function (4-22) would increase. In the case when the value of the function (g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2 is greater than the value of the function r.sup.2 /(p.sup.2 +q.sup.2 (.parallel..gradient.u.sup..sigma. (x,y,t).parallel..sup.2 +.parallel..gradient.v.sup..sigma. (x,y,t).parallel..sup.2)), the fact that the summand of the quadratic form (4-21) associated with the optical flow constraint g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.) is negative implies that the function (4-22) would decrease. Therefore, the function (4-22) achieves its maximal value ##EQU32## when the value of the function (g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2 is equal to the value of the function r.sup.2 /(p.sup.2 +q.sup.2 (.parallel..gradient.u.sup..sigma. (x,y,t).parallel..sup.2 +.parallel..gradient.v.sup..sigma. (x,y,t).parallel..sup.2)).

It may also happen that the values of the functions .parallel..gradient.u.sup..sigma. (x,y,t).parallel. and .parallel..gradient.v.sup..sigma. (x,y,t).parallel. become large. Under these conditions the following values are reduced: the maximal value of the function (4-22) specified by the relation (4-24), the value r.sup.2 /(p.sup.2 +q.sup.2 (.parallel..gradient.u.sup..sigma. (x,y,t).parallel..sup.2 +.parallel..gradient.v.sup..sigma. (x,y,t).parallel..sup.2)) of the function (g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2 where this maximal value is achieved. The above situation may occur in the following cases.

1. The occluding object may artificially constrain the estimate (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) in the direction orthogonal to the occluding boundary. Then the optical flow constraint may react by shifting the estimate (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) of the velocity vector in the direction parallel to the occluding boundary to reduce the absolute value of the function g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.). This, in turn, may increase the values of the functions .parallel..gradient.u.sup..sigma. (x,y,t).parallel. and .parallel..gradient.v.sup..sigma. (x,y,t).parallel..

2. The estimate (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) may diverge as the result of nonlinearity of the optical flow constraint or as the result of the noise present in the images.

The directional smoothness constraints (s,.gradient.u.sup..sigma. (x,y,t)) and (s,.gradient.v.sup..sigma. (x,y,t)) in the direction s.epsilon.S are represented in the system of nonlinear equations (4-19) by the expression ##EQU33## With respect to this expression the following cases can be considered:

1. The image point (x,y) does not cross the occluding boundary in the direction s. Then the value of the function ((s,.gradient.u.sup..sigma. (x,y,t)).sup.2 +(s,.gradient.v.sup..sigma. (x,y,t)).sup.2) may become small. In this case the value of the function (4-25) would approach the value of the function ##EQU34## This, in turn, would make each of the two summands of the system of nonlinear equations (4-19), one associated with the directional smoothness constraint (s,.gradient.u.sup..sigma. (x,y,t)), and the other associated with the directional smoothness constraint (s,.gradient.v.sup..sigma. (x,y,t)), behave as a linear function of the corresponding constraint with proportionality coefficient equal to the constant value 1/a.sup.2.

2. The image point (x,y) approaches the occluding boundary in the direction s. Then the value of the function ((s,.gradient.u.sup..sigma. (x,y,t)).sup.2 +(s,.gradient.v.sup..sigma. (x,y,t)).sup.2) may increase. In the case when the value of the function ((s,.gradient.u.sup..sigma. (x,y,t)).sup.2 +(s,.gradient.v.sup..sigma. (x,y,t)).sup.2) is less than the value of the function a.sup.2 /(c.sup.2 +b.sup.2 (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2), the fact that the summand of the quadratic form (4-21) associated with the directional smoothness constraints (s,.gradient.u.sup..sigma. (x,y,t)) and (s,.gradient.v.sup..sigma. (x,y,t)) is positive implies that the function (4-25) would increase. In the case when the value of the function ((s,.gradient.u.sup..sigma. (x,y,t)).sup.2 +(s,.gradient.v.sup..sigma. (x,y,t)).sup.2) is greater than the value of the function a.sup.2 /(c.sup.2 +b.sup.2 (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2), the fact that the summand of the quadratic form (4-21) associated with the directional smoothness constraints (s,.gradient.u.sup..sigma. (x,y,t)) and (s,.gradient.v.sup..sigma. (x,y,t)) is negative implies that the function (4-25) would decrease. Therefore, the function (4-25 ) achieves its maximal value ##EQU35## when the value of the function ((s,.gradient.u.sup..sigma. (x,y,t)).sup.2 +(s,.gradient.v.sup..sigma. (x,y,t)).sup.2) is equal to the value of the function a.sup.2 /(c.sup.2 +b.sup.2 (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2).

In the case when the values of some of the functions (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)), g.sub.t .epsilon.G.sub.t are large the following values are reduced: the maximal value of the function (4-25) specified by the relation (4-27), the value a.sup.2 /(c.sup.2 +b.sup.2 (s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2) of the function ((s,.gradient.u.sup..sigma. (x,y,t)).sup.2 +(s,.gradient.v.sup..sigma. (x,y,t)).sup.2) where this maximal value is achieved. Such an action reflects the fact that in this case the image point (x,y) is likely to cross the occluding boundary in the direction s if the value of the function ((s,.gradient.u.sup..sigma. (x,y,t)).sup.2 +(s,.gradient.v.sup..sigma. (x,y,t)).sup.2) specifying the degree of the smoothness of the estimate of the velocity vector field in such a direction s is large.

As has been stated earlier the role of the regularization constraint is to discriminate between the parametric vector fields (u,v) giving the same optimal values to the weighted average of the optical flow and the directional smoothness constraints without causing any significant changes in these values. This is achieved by assigning a small value to the parameter .gamma. appearing in the system of nonlinear equations (4-19). For the sake of simplicity of the analysis given below, we shall ignore the regularization constraint by assuming that the parameter .gamma. is equal to zero while keeping in mind that the solution of the system of nonlinear equations (4-19) is locally unique.

Taking the above into account, the process of solving the system of nonlinear equations (4-19) described in the next section can be regarded as an attempt to achieve an equal balance between the functions representing two antagonistic sets of constraints: optical flow constraints and directional smoothness constraints. This, in part, is the consequence of the fact that the quadratic form arising from the linearization (5-16) of the system of nonlinear equations (4-19) described in the next section is positive definite.

It is not difficult to see that by increasing the values of the parameters r.sup.2,p.sup.2,q.sup.2 appearing in the relation (4-12), we decrease the values of the weight functions .alpha..sub.g.sbsb.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,.parallel..gradient.u.sup..sigma. .parallel.,.parallel..gradient.v.sup..sigma. .parallel.), g.sub.t .epsilon.G.sub.t. This would shift the balance in the system of nonlinear equations (4-19) toward the directional smoothness constraints. To restore the equal balance between the optical flow and the directional smoothness constraints the system of nonlinear equations (4-19) would react in such a way that the absolute values of the functions (s,.gradient.u.sup..sigma. (x,y,t)), (s,.gradient.v.sup..sigma. (x,y,t)), s.epsilon.S corresponding to the directional smoothness constraints are reduced, while the absolute values of the functions g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.), g.sub.t .epsilon.G.sub.t corresponding to the optical flow constraints are increased. Likewise, by increasing the values of the parameters a.sup.2,c.sup.2,b.sub.g.sbsb.t.sup.2,g.sub.t .epsilon.G.sub.t appearing in the relation (4-13), we decrease the values of the weight functions .beta..sub.s.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.,(s,.gradient.u.sup..sigma.),(s,.gradient.v.sup..sigma.)), s.epsilon.S. This, in turn, would shift the balance in the system of nonlinear equations (4-19) toward the optical flow constraints. To restore equal balance between the optical flow and the directional smoothness constraints the system of nonlinear equations (4-19) would react in the way that would reduce the absolute values of the functions g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.), g.sub.t .epsilon.G.sub.t corresponding to the optical flow constraints and would increase the absolute values of the functions (s,.gradient.u.sup..sigma. (x,y,t)), (s,.gradient.v.sup..sigma. (x,y,t)), s.epsilon.S corresponding to the directional smoothness constraints.

Taking into account the earlier discussion related to the behavior of the optical flow and the directional smoothness constraints, we can conclude the following. The parameters r.sup.2 and a.sup.2 have their greatest impact on the solution of the system of nonlinear equations (4-19) at the image points that are away from the occluding boundary, while the parameters p.sup.2,q.sup.2,c.sup.2,b.sub.g.sbsb.t.sup.2, g.sub.t .epsilon.G.sub.t mostly influence the solution at the image points that are near the occluding boundary. The parameter r.sup.2 defines the coefficient of proportionality for the functions g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.), g.sub.t .epsilon.G.sub.t specifying the optical flow constraints, while the parameter a.sup.2 determines the proportionality coefficient for the functions (s,.gradient.u.sup..sigma. (x,y,t)), (s,.gradient.v.sup..sigma. (x,y,t)), s.epsilon.S specifying the directional smoothness constraints. The combination (4-24) of the parameters r.sup.2, p.sup.2, and q.sup.2 determines the upper limit for the influence of the optical flow constraints on the solution of the system of nonlinear equations (4-19), while the combination (4-27) of the parameters a.sup.2, c.sup.2, and b.sub.g.sbsb.t.sup.2, g.sub.t .epsilon.G.sub.t determines the upper limit for the influence of the directional smoothness constraints on the solution of the system of nonlinear equations (4-19).

Earlier the estimate {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} was defined as a solution of the system of nonlinear equations (4-19), which can be expressed in the following form: ##EQU36## Here, for the sake of compactness of the notation, the arguments of the functions appearing in the system of nonlinear equations (4-19) have been omitted. Also, the symbols g.sub.t.sup..sigma., g.sub.tu.sup..sigma., g.sub.tv.sup..sigma., g.sub.t.sup..sigma. have been used instead of the symbols g.sub.t.sup..sigma., g.sub.tu.sup..sigma., g.sub.tv.sup..sigma., g.sub.t.sup..sigma. to indicate that at every image point (x,y).epsilon..OMEGA. the corresponding functions depend on the estimate (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) of the velocity vector (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).

To compute the estimate of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.,.sigma..epsilon.[1,.infin.)} a restriction is imposed ourselves to a finite subset .sigma..sub.1, . . . , .sigma..sub.n of the parameter values from the interval [1,.infin.). Let .sigma..sub.0 be the parameter with the value equal to .infin.. Suppose that the value of the parameter .sigma..sub.n is equal to 1; and that the sequence of the parameters .sigma..sub.0, .sigma..sub.1, . . . , .sigma..sub.n is not too sparse, and is decreasing in values. Then the estimate of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma.=.sigma..sub.1, . . . , .sigma..sub.n } is obtained by successively solving the system of nonlinear equations (5-1) for each value of the parameter of .sigma.=.sigma..sub.1, . . . , .sigma..sub.n, starting with the parameter value .sigma.=.sigma..sub.0 where the solution is identically equal to a given vector constant (u.sup..infin.,v.sup..infin.), and proceeding in the direction decreasing with respect to the parameter values using the solution of the system of nonlinear equations (5-1) obtained for the preceding parameter value as the initial estimate.

Let .delta.(.sigma.) be a function of the parameter .sigma. defined on the set of the parameter values .sigma.=.sigma..sub.1, . . . , .sigma..sub.n by the following relation:

.delta.(.sigma..sub.i)=.sigma..sub.i-1, i=1, . . . , n. (5-2) If the weight constants .gamma..sup..sigma., .sigma.=.sigma..sub.1, . . . , .sigma..sub.n are defined as in the relation .gamma.=.gamma..sup..sigma. (.delta.(.sigma.)-.sigma.), (5-3) then the regularization constraints ##EQU37## can, respectively, be approximated as ##EQU38## for every image point (x,y).epsilon..OMEGA., and for every parameter value .sigma.=.sigma..sub.1, . . . , .sigma..sub.n. Under these conditions the system of nonlinear equations (5-1) of the unknown estimate of the parametric velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA., .sigma.=.sigma..sub.1, . . . , .sigma..sub.n } takes the form ##EQU39## Then the computational process for estimating the parametric velocity vector field can be more precisely described as follows: for every parameter value .sigma.=.sigma..sub.1, . . . , .sigma..sub.n, using the previously obtained estimate {(u.sup..delta.(.sigma.) (x,y,t),v.sup..delta.(.sigma.) (x,y,t)).vertline.(x,y).epsilon..OMEGA.} of the velocity vector field corresponding to the parameter value .delta.(.sigma.) as the initial estimate, find the estimate {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.} of the velocity vector field corresponding to the parameter value .sigma. that is the closest to the initial estimate by solving the system of nonlinear equations (5-6).

For a fixed value of the parameter .sigma.=.sigma..sub.1, . . . , .sigma..sub.n, the system of nonlinear equations (5-6) relative to the unknown estimate (u.sup..sigma.,v.sup..sigma.).tbd.{(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.} of the velocity vector field corresponding to the parameter value .sigma., defines, for every image point (x,y).epsilon..OMEGA., a set of relations between the components of the velocity vector (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) and its partial derivatives, which can be expressed in the form

F.sup..sigma. (u.sup..sigma.,v.sup..sigma.)=0. (5-7) To find the estimate (u.sup..sigma.,v.sup..sigma.) of the velocity vector field corresponding to the parameter value .sigma. as an estimate of the solution of the system of nonlinear equations (5-6) corresponding to the parameter value .sigma. that is the closest to the initial estimate (u.sup..delta.(.sigma.),v.sup..delta.(.sigma.)) the iterative updating scheme can be used. In this scheme, starting with the initial estimate (u.sup..delta.(.sigma.),v.sup..delta.(.sigma.)), which is taken as the current estimate (u.sup..sigma.,v.sup..sigma.) of the solution of the system of nonlinear equations (5-6) corresponding to the parameter value .sigma., the improved estimate (u.sub.+.sup..sigma.,v.sub.+.sup.94 ).tbd.{(u.sub.+.sup..sigma. (x,y,t),v.sub.+.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.} of the solution of the system of nonlinear equations (5-6) corresponding to the parameter value .sigma. is determined by the step (u.sub.+.sup..sigma.,v.sub.+.sup..sigma.)=(u.sup..sigma.,v.sup..sigma.)+.omega.(.DELTA.u.sup..sigma., .DELTA.v.sup..sigma.), (5-8) where the scalar parameter .omega. defines the length of the step, while the vector field (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).tbd.{(.DELTA.u.sup..sigma. (x,y,t),.DELTA.v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.} defines the direction of the step. The step length .omega..epsilon.(0,1] is selected in such a way that the following function is minimized: .function..sup..sigma. (.omega.)=F.sup..sigma. (u.sup..sigma. +.omega..DELTA.u.sup..sigma.,v.sup..sigma. +.omega..DELTA.v.sup..sigma.).sup.T F.sup..sigma. (u.sup..sigma. +.omega..DELTA.u.sup..sigma.,v.sup..sigma. +.omega..DELTA.v.sup..sigma.). (5-9) The step direction (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) is selected in such a way that the function (5-9) becomes steadily decreasing for sufficiently small values of .omega..epsilon.(0,1]. The improved estimate (u.sub.+.sup..sigma.,v.sub.+.sup..sigma.) is taken as the current estimate (u.sup..sigma.,v.sup..sigma.) of the solution of the system of nonlinear equations (5-6) corresponding to the parameter value .sigma., and the step (5-8) of the iterative updating scheme is repeated to determine the next improved estimate (u.sub.+.sup..sigma.,v.sub.+.sup..sigma.) of the solution of the system of nonlinear equations (5-6) corresponding to the parameter value .sigma.. This process is continued until the appropriate criteria are met, in which case the improved estimate (u.sub.+.sup.94 ,v.sub.+.sup..sigma.) is taken as the estimate (u.sup..sigma.,v.sup..sigma.) of the velocity vector field corresponding to the parameter value .sigma..

If the vector field (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) is not too large, then in the neighborhood of the vector field (u.sup..sigma.,v.sup..sigma.) the nonlinear operator F.sup..sigma. can be linearly expanded as

F.sup..sigma. (u.sup..sigma. +.DELTA.u.sup..sigma.,v.sup..sigma. +.DELTA.v.sup..sigma.).apprxeq.F.sup..sigma. (u.sup..sigma.,v.sup..sigma.)+J.sup..sigma. (u.sup..sigma.,v.sup..sigma.,.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.). (5-10) Here the linear and bounded with respect to the vector field (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) operator J.sup..sigma. (u.sup..sigma.,v.sup..sigma.,.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) is the Jacobian of the nonlinear operator F.sup..sigma. (u.sup..sigma.,v.sup..sigma.). The operator J.sup..sigma. (u.sup..sigma.,v.sup..sigma.,.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) can be more conveniently expressed in the form J.sup..sigma. (u.sup..sigma.,v.sup..sigma.,.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).tbd.J.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.). (5-11)

In the Newton method [22] the vector field (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) is defined as a solution of the system of linear equations

J.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)=-F.sup..sigma. (u.sup..sigma.,v.sup..sigma.). (5-12) Then the improved estimate (u.sub.+.sup..sigma.,v.sub.+.sup..sigma.) of the solution of the system of nonlinear equations (5-6) is defined as in the relation (5-8) where the scalar parameter .omega. is taken to be equal to 1.

The improved estimate (u.sub.+.sup..sigma.,v.sub.+.sup..sigma.) of the solution of the system of nonlinear equations (5-6), which is obtained as the result of solving the system of linear equations (5-12) and then applying the relation (5-8), is not necessarily a better estimate. The reason behind it comes from the fact that the Jacobian J.sup..sigma. (u.sup..sigma.,v.sup..sigma.) of the nonlinear operator F.sup..sigma. (u.sup..sigma.,v.sup..sigma.) is nonsymmetric and ill-conditioned, so that the system of linear equations (5-12) cannot be reliably solved for the vector field (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.). The last is the consequence of the fact that the regularization constraints can only be effective in resolving the ambiguities of the system of nonlinear equations (5-6) at the image points that are beyond the occluding boundaries where the quadratic form appearing in the relation (4-21) is positive definite. The regularization constraints become ineffective in resolving the ambiguities of the system of nonlinear equations (5-6) at the image points that are near the occluding boundaries where the quadratic form appearing in the relation (4-21) is not positive definite. This observation suggests that an alternative method should be employed.

In the quasi-Newton method of solving the system of nonlinear equations (5-6) the system of linear equations (5-12) is replaced with the system of linear equations

M.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)=-F.sup..sigma. (u.sup..sigma.,v.sup..sigma.), (5-13) which combines the globally convergent strategy of the well-conditioned method with the fast local strategy of the Newton method in a way that derives the benefits of both. This quasi-Newton method defines the regularization of the system of nonlinear equations (5-6) that is meaningful for the problem of estimating the parametric velocity vector field. The linear and bounded with respect to the vector field (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) operator M.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) is the approximation to the Jacobian J.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) of the nonlinear operator F.sup..sigma. (u.sup..sigma.,v.sup..sigma.) defined as follows.

Let F.sub..theta..sup..sigma. (u.sup..sigma.,v.sup..sigma.,u.sup..sigma.,v.sup..sigma.), .theta..epsilon.[0,1], be a family of nonlinear operators of the invariable copy (u.sup..sigma.,v.sup..sigma.) and of the variable copy (u.sup..sigma.,v.sup..sigma.) of the estimate of the velocity vector field corresponding to the parameter value .sigma. defined by the relation ##EQU40## where the symbols g.sub.t.sup..sigma., g.sub.tu.sup..sigma., g.sub.tv.sup..sigma., g.sub.t.sup..sigma. are used to indicate that at every image point (x,y).epsilon..OMEGA. the corresponding functions depend on the invariable copy (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) of the estimate of the velocity vector (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)), while the symbol g.sub.t.sup..sigma. is used to indicate that at every image point (x,y).epsilon..OMEGA. the corresponding function depends on the variable copy (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)) of the estimate of the velocity vector (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)). If the variable copy (u.sup..sigma.,v.sup..sigma.) of the estimate of the velocity vector field is identically equal to the invariable copy (u.sup..sigma.,v.sup..sigma.) of the estimate of the velocity vector field, then for every .theta..epsilon.[0,1] the nonlinear operator F.sub..theta..sup..sigma. (u.sup..sigma.,v.sup..sigma.,u.sup..sigma.,v.sup..sigma.) is identically equal to the nonlinear operator F.sup..sigma. (u.sup..sigma.,v.sup..sigma.). The parameter .theta. defines the degree of the feedback relaxation of the optical flow and the directional smoothness constraints through the variable copy (u.sup..sigma.,v.sup..sigma.) of the estimate of the velocity vector field.

Let M.sub..theta..sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.), .theta..epsilon.[0,1] be a family of the linear and bounded with respect to the vector field (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) operators, where for each .theta..epsilon.[0,1] the operator M.sub..theta..sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) is the Jacobian of the nonlinear operator F.sub..theta..sup..sigma. (u.sup..sigma.,v.sup..sigma.,u.sup..sigma.,v.sup..sigma.), considered as the function of the vector field (u.sup..sigma.,v.sup..sigma.), and depending on the vector field (u.sup..sigma.,v.sup..sigma.) as on the parameter, under the conditions that the vector field (u.sup..sigma.,v.sup..sigma.) is identically equal to the vector field (u.sup..sigma.,v.sup..sigma.). Then the operator M.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) is defined as the member of the family M.sub..theta..sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.), .theta..epsilon.[0,1] satisfying the following requirements:

1. As a part of the globally convergent strategy, the vector field (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) obtained as the solution of the system of linear equations (5-13) is required to be a descent direction for the functional h(u.sup..sigma.,v.sup..sigma.,u.sup..sigma.,v.sup..sigma.), considered as the function of the vector field (u.sup..sigma.,v.sup..sigma.), and depending on the vector field (u.sup..sigma.,v.sup..sigma.) as on the parameter, under the conditions that the vector field (v.sup..sigma.,v.sup..sigma.) be identically equal to the vector field (u.sup..sigma.,v.sup..sigma.). In other words, the first-order partial with respect to the variable .omega. derivative of the functional h(u.sup..sigma.,v.sup..sigma.,u.sup..sigma. +.omega..DELTA.u.sup..sigma.,v.sup..sigma. +.omega..DELTA.v.sup..sigma.) considered as the function of the scalar variable .omega. is required to be negative when the value of such variable .omega. is equal to 0. It is not difficult to see that the value of this derivative is equal to

-F.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.T M.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-1 F.sup..sigma. (u.sup..sigma.,v.sup..sigma.). (5-15) The requirement for the functional (5-15) to be negative can be satisfied by in turn requiring the operator M.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) to be positive definite.

2. As a part of the fast local strategy, the operator M.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) is required to be the nearest possible approximation to the Jacobian J.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma., .DELTA.v.sup..sigma.) of the nonlinear operator F.sup..sigma. (u.sup..sigma.,v.sup..sigma.).

The operator M.sub..theta..sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) from the family corresponding to the parameter value .theta.=0.5 satisfies the above requirements and is taken as the operator M.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) appearing in the system of linear equations (5-13). With this choice of the operator M.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) the system of linear equations (5-13) has the following explicit form ##EQU41##

Next, the finite-difference discretizations of the system of nonlinear equations (5-6) as well as the system of linear equations (5-16) are considered. The method of constructing the estimate of the solution of the system of nonlinear equations (5-6) is extended to the finite-difference discretization of this system.

Although the theoretical analysis for the computation of the estimate of the parametric velocity vector field can be effectively performed in the case when the combined generalized initial image functions .GAMMA..sub.g, g.epsilon.G are defined over the entire set of testing functions .PHI.(R.sup.3), while the estimate (u(x,y,t),v(x,y,t)) of the velocity vector (u,(x,y,t),v(x,y,t)) is defined for every image point (x,y).epsilon..OMEGA. and for every time moment t.epsilon.T, the actual computation can only be accomplished on finite subsets of the above data sets. In this section a finite-difference discretization of the system of nonlinear equations (5-6) is described. The estimate of the parametric velocity vector field is then determined as a solution of the resulted system of a finite number of nonlinear equations.

Given a generalized function F defined on the locally convex linear topological space .PHI.(R.sup.3), and given a testing function .phi..epsilon..PHI.(R.sup.3) the "convolution" of the generalized function F and the testing function .phi. is the infinitely differentiable function (F*.phi.) defined on the space R.sup.3 by the relation

(F*.phi.)(x,y,t)=F(.phi..sub.x,y,t), .phi..sub.x,y,t (x,y,t).tbd..phi.(x-x,y-y,t-t), (x,y,t).epsilon.R.sup.3, (6-1) for every (x,y,t).epsilon.R.sup.3.

Let F be a generalized function, and let .chi..epsilon..PHI.(R.sup.3) be a given fixed non-negative testing function such that ##EQU42## For instance, take ##EQU43## where .psi. is the function described later. The testing function .chi. will be called the "sampling function". Then the "regularization" of the generalized function F through the sampling function .chi. is defined on the space R.sup.3 as the infinitely differentiable function (F*.chi.), which is obtained as the convolution of the generalized function F and the sampling function .chi..

Let {t'.sub.k }.sub.k=0.sup.K' be an increasing sequence of time moments containing the sequence of time moments {t".sub.k }.sub.k=0.sup.K". Assume that for every t.epsilon.{t'.sub.k }.sub.k=0.sup.K' the regularizations .gamma..sub.86 .tbd.(.GAMMA..sub..xi. *.chi.) of the generalized functions .GAMMA..sub..xi. associated with the irradiance image functions .xi., .xi..epsilon..XI. through the sampling function .chi. and the regularizations .gamma..sub..eta. .tbd.(.GAMMA..sub..eta. *.chi.) of the generalized functions .GAMMA..sub..eta. associated with the feature image functions .eta., .eta..epsilon.H through the sampling function .chi. are given at the points (x,y) belonging to the following subset of the image plane:

R.sup.2 (h'.sub..sigma.,1 h'.sub..sigma.,2)={(x',y').vertline.(x',y')=i.sub.1 h'.sub..sigma.,1 +i.sub.2 h'.sub..sigma.,2,i.sub.1,i.sub.2 .epsilon.Z}. (6-2) Here Z is the set of integers; h'.sub..sigma.,1, h'.sub..sigma.,2 are two dimensional real vectors; and i.sub.1 h'.sub..sigma.,1 +i.sub.2 h'.sub..sigma.,2 is a linear combination of such vectors with integer coefficients i.sub.1, i.sub.2. The points from the subset R.sup.2 (h'.sub..sigma.,1,h'.sub..sigma.,2) will be called the "sampling points" corresponding to the parameter .sigma.; the functions .gamma..sub..xi., .xi..epsilon..XI. will be called the "sampled irradiance image functions"; and the functions .gamma..sub..eta., .eta..epsilon.H will be called the "sampled feature image functions".

The generalized functions .GAMMA..sub..xi..sigma. associated with the sampled irradiance image functions .gamma..sub..xi., .xi..epsilon..XI. and the generalized functions .GAMMA..sub..eta.,.sigma. associated with the sampled feature image functions .gamma..sub..eta., .eta..epsilon.H are defined on the set of testing functions .PHI.(R.sup.3) as follows. The value of the generalized function .GAMMA..sub..xi.,.sigma. associated with the sampled irradiance image function .gamma..sub..xi. at the testing function .phi..epsilon..PHI.(R.sup.3) is given by the relation ##EQU44## The value of the generalized function .GAMMA..sub..eta.,.sigma. associated with the sampled feature image function .gamma..sub..eta. at the testing function .phi..epsilon..PHI.(R.sup.3) is given be the relation ##EQU45##

Let .lambda..tbd.{.lambda..sub..xi.,.lambda..sub.72 .vertline..xi..epsilon..XI., .eta..epsilon.H} be a set of real-valued constants; let g(.lambda.) be an index attached to the set .lambda.; and let m.sub.g(.lambda.) .tbd.{m.sub..xi.,x,m.sub..xi.,y,m.sub..xi.,t, m.sub..eta., x,m.sub..eta.,y,m.sub..eta.,t .vertline..xi..epsilon..XI., .eta..epsilon.H} be a set of non-negative integer constants corresponding to the set of constants .lambda.. Then the "combined generalized sampled image function" corresponding to the index g.tbd.g(.lambda.) is the generalized function .GAMMA..sub.g,.sigma. defined on the locally convex linear topological space .PHI.(R.sup.3). The value of the combined generalized sampled image function .GAMMA..sub.g,.sigma. at the testing function .phi..epsilon..PHI.(R.sup.3) is given by the relation ##EQU46## and is called the "observation" of the combined generalized sampled image function .GAMMA..sub.g,.lambda. on the testing function .phi..

Assume the restriction to the case when each subset R.sup.2 (h'.sub..sigma.,1,h'.sub..sigma.,2), .sigma.=.sigma..sub.0, .sigma..sub.1, . . . , .sigma..sub.n is a square grid of points on the image plane. In other words, it is a requirement that the vectors h'.sub..sigma.,1, h'.sub..sigma.,2 are of equal length and are orthogonal to each other for every .sigma.=.sigma..sub.0, .sigma..sub.1, . . . , .sigma..sub.n. Also, it is a requirement that such grids of points constitute a pyramid on the image plane which means that the following relation is satisfied:

R.sup.2 (h'.sub..sigma..sbsb.0.sub.,1,h'.sub..sigma..sbsb.0.sub.,2).OR right.R.sup.2 (h'.sub..sigma..sbsb.1.sub.,1,h'.sub..sigma..sbsb.1.sub.,2).OR right. . . . .OR right.R.sup.2 (h'.sub..sigma..sbsb.n.sub.,1,h'.sub..sigma..sbsb.n.sub.,2). (6-6)

The "grid rotation method" is the sampling method that complies with the above requirements. It is commonly used to obtain numerical solutions of partial differential equations. In this method the vectors h'.sub..sigma..sbsb.0.sub.,1, h'.sub..sigma..sbsb.0.sub.,2 are selected to be of equal length and orthogonal to each other; then for every k=1, . . . , n the vectors h'.sub..sigma..sbsb.k.sub.,1, h'.sub..sigma..sbsb.k.sub.,2 are defined as in the relations

h'.sub..sigma..sbsb.k.sub.,1 =0.5(h'.sub..sigma..sbsb.k-1.sub.,1 +h'.sub..sigma..sbsb.k-1.sub.,2), (6-7) h'.sub..sigma..sbsb.k.sub.,2 =0.5(h'.sub..sigma..sbsb.k-1.sub.,1 -h'.sub..sigma..sbsb.k-1.sub.,2), (6-8) It is not difficult to see that the above defined vectors h'.sub..sigma..sbsb.k.sub.,1, h'.sub..sigma..sbsb.k.sub.,2 are of equal length and are orthogonal to each other for every k=0, 1, . . . , n, and that the relation 6-6) is satisfied.

Let .kappa.' be a positive integer constant and .psi..epsilon..PHI.(R.sup.3) be a given measurement function. Then for every time moment t.epsilon.{t".sub.k }.sub.k=0.sup.K", and for every parameter value .sigma.=.sigma..sub.0, .sigma..sub.1, . . . , .sigma..sub.n the image function g.sup..sigma. (x,y,t) and its spatial partial derivatives g.sub.x.sup..sigma. (x,y,t), g.sub.y.sup..sigma. (x,y,t) are defined on the set .OMEGA. as follows.

First, assume that the image point (x,y) belongs to the subset .OMEGA."(h".sub..sigma.,1,h".sub..sigma.,2).andgate.R.sup.2 (h".sub..sigma.,1,h".sub..sigma.,2). Here the two-dimensional real vectors h".sub..sigma.,1, h".sub..sigma.,2 are defined as in the relations

h".sub..sigma.,1 =.kappa.'h'.sub..sigma.,1, (6-9) h".sub..sigma.,2 =.kappa.'h'.sub..sigma.,2, (6-10) while the subset .OMEGA."(h".sub..sigma.,1,h".sub..sigma.,2).OR right.R.sup.2 is defined as in the relation ##EQU47## Then for every g.epsilon.G the value of the component g.sup..sigma. (x,y,t) of the image function g.sup..sigma. (x,y,t) is determined as the observation .GAMMA..sub.g,.sigma. (.psi..sub.x,y,t.sup..sigma.) of the combined generalized sampled image function .GAMMA..sub.g,.sigma. on the testing function .psi..sub.x,y,t.sup..sigma. defined by the relation (2-4); and the values of the components g.sub.x.sup..sigma. (x,y,t), g.sub.y.sup..sigma. (x,y,t) of the partial derivatives g.sub.x.sup..sigma. (x,y,t), g.sub.y.sup..sigma. (x,y,t) of the image function g.sup..sigma. (x,y,t) are given as the observations ##EQU48## of the combined generalized sampled image function .GAMMA..sub.g,.sigma. on the defined by the relations (2-5), (2-6) testing functions ##EQU49##

Let (x,y) be any image point from the set .OMEGA., then for some integer numbers i.sub.1, i.sub.2 and for some real number .theta..sub.1, .theta..sub.2 from the interval [0,1] the following relation is satisfied:

(x,y)=(i.sub.1 +.theta..sub.1)h".sub..sigma.,1 +(i.sub.2 +.theta..sub.2)h".sub..sigma.,2. (6-12) Then the values of the image function g.sup..sigma. (x,y,t) and its partial derivatives g.sub.x.sup..sigma. (x,y,t), g.sub.y.sup..sigma. (x,y,t) are obtained through the bilinear interpolations g.sup..sigma. (x,y,t)=(1-.theta..sub.2)((1-.theta..sub.1)g.sup..sigma. (x.sub.1,1,y.sub.1,1,t)+.theta..sub.1 g.sup..sigma. (x.sub.1,2,y.sub.1,2,t))+.theta..sub.2 ((1-.theta..sub.1)g.sup..sigma. (x.sub.2,1,y.sub.2,1,t)+.theta..sub.1 g.sup..sigma. (x.sub.2,2,y.sub.2,2,t)), (6-13) g.sub.x.sup..sigma. (x,y,t)=(1-.theta..sub.2)((1-.theta..sub.1)g.sub.x.sup..sigma. (x.sub.1,1,y.sub.1,1,t)+.theta..sub.1 g.sub.x.sup..sigma. (x.sub.1,2,y.sub.1,2,t))+.theta..sub.2 ((1-.theta..sub.1)g.sub.x.sup..sigma. (x.sub.2,1,y.sub.2,1,t)+.theta..sub.1 g.sub.x.sup..sigma. (x.sub.2,2,y.sub.2,2,t)), (6-14) g.sub.y.sup..sigma. (x,y,t)=(1-.theta..sub.2)((1-.theta..sub.1)g.sub.y.sup..sigma. (x.sub.1,1,y.sub.1,1,t)+.theta..sub.1 g.sub.y.sup..sigma. (x.sub.1,2,y.sub.1,2,t))+.theta..sub.2 ((1-.theta..sub.1)g.sub.y.sup..sigma. (x.sub.2,1,y.sub.2,1,t)+.theta..sub.1 g.sub.y.sup..sigma. (x.sub.2,2,y.sub.2,2,t)), (6-15) of the values of the image function g.sup..sigma. (x,y,t) and, respectively, its partial derivatives g.sub.x.sup..sigma. (x,y,t), g.sub.y.sup..sigma. (x,y,t) at the points from the subset .OMEGA."(h".sub..sigma.,1,h".sub..sigma.,2).andgate.R.sup.2 (h".sub..sigma.,1,h".sub..sigma.,2) given as in the relations (x.sub.1,1,y.sub.1,1)=i.sub.1 h".sub..sigma.,1 +i.sub.2 h".sub..sigma.,2, (6-16) (x.sub.1,2,y.sub.1,2)=(i.sub.1 +1)h".sub..sigma.,1 +i.sub.2 h".sub..sigma.,2, (6-17) (x.sub.2,1,y.sub.2,1)=i.sub.1 h".sub..sigma.,1 +(i.sub.2 +1)h".sub..sigma.,2, (6-18) (x.sub.2,2,y.sub.2,2)=(i.sub.1 +1)h".sub..sigma.,1 +(i.sub.2 +1)h".sub..sigma.,2. (6-19)

Let .kappa." be a positive integer constant, and let h.sub..sigma.,1, h.sub..sigma.,2 be two-dimensional real vectors defined as in the relations

h.sub..sigma.,1 =.kappa."h".sub..sigma.,1, (6-20) h.sub..sigma.,2 =.kappa."h".sub..sigma.,2, (6-21) while the subset .OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2).OR right..OMEGA. be defined as in the relation .OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2)=.OMEGA..andgate.R.sup.2 (h.sub..sigma.,1,h.sub..sigma.,2). (6-22) The finite-difference discretizations of the system of nonlinear equations (5-6) and of the system of linear equations (5-16), both corresponding to the parameter value .sigma.=.sigma..sub.0, .sigma..sub.1, . . . , .sigma..sub.n, are defined on the set .OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2) as, respectively, the system of nonlinear equations relative to the unknown estimate of the velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2)} corresponding to the parameter value .sigma. and its finite differences, and the system of linear equations relative to the unknown vector field {(.DELTA.u.sup..sigma. (x,y,t),.DELTA.v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2)} corresponding to the parameter value .sigma. and its finite differences as follows.

The expressions appearing under the integral signs over the set S in the system on nonlinear equations (5-6) and in the system of linear equations (5-16) have the general form

(s,.gradient..function.(x,y,t,s)), (6-23) where .function.(x,y,t,s) is a function of the variables (x,y).epsilon..OMEGA., t.epsilon.T, s.epsilon.S, which is continuously differentiable with respect to the variables (x,y).epsilon..OMEGA. for every fixed value of the variables t.epsilon.T, s.epsilon.S.

The unit circle making up the set S is replaced by a finite subset, denoted by the same symbol S, of the set R.sup.2 (h.sub..sigma.,1,h.sub..sigma.,2) having the following properties: the set S does not contain the origin, and for every element s belonging to the set S the element -s also belongs to the set S. The measure ds is replaced with the point measure associating the value 0.5 to every element s.epsilon.S.

As a consequence of the above actions an integral over the set S of values of any expression of the form (6-23) becomes a finite sum over the set S of the values of that expression; in other words, the following relation is satisfied: ##EQU50##

Given a function .function.(x,y,t,s) that is continuously differentiable with respect to the variables (x,y).epsilon..OMEGA. for every t.epsilon.T, s.epsilon.S, and given an element s.epsilon.S, the expression (6-23) is defined on the set R.sup.2 (h.sub..sigma.,1,h.sub..sigma., 2) in terms of the finite-difference approximations to the partial derivatives as

(s,.gradient..function.(x,y,t,s))=.rho..sub.s (.function.(x+0.5s.sub.x,y+0.5s.sub.y,t,s)-.function.(x-0.5s.sub.x,y-0.5s.sub.y,t,s)), (6-25) where .rho..sub.s is a real positive constant depending on the length of the vector s as a parameter, and the values s.sub.x, s.sub.y are the components of the vector s.tbd.(s.sub.x,s.sub.y).

Taking into account the relation (6-25) and the fact that the set S has the property that for every element s.epsilon.S the element -s.epsilon.S, then, in the case when .function.(x,y,t,s) is an odd function with respect to the variable s, the sum appearing in the relation (6-24) can be expressed as ##EQU51##

It is assumed that the set G.sub.t is finite and that dg.sub.t is the point measure associating the value .rho.(g.sub.t) to the element g.sub.t .epsilon.G.sub.t. Here .rho.(g.sub.t) is a real positive constant for every g.sub.t .epsilon.G.sub.t.

Based on the above discussion, the system of nonlinear equations (5-6) can be defined on the set R.sup.2 (h.sub..sigma.,1,h.sub..sigma.,2) as ##EQU52## while the system of linear equations (5-16) can be defined on the set R.sup.2 (h.sub..sigma.,1,h.sub..sigma.,2) as ##EQU53## The arguments of the functions appearing in the relations (6-27), (6-28) have been omitted for the sake of simplicity of notation. The functions appearing under the summation over the set S sign are evaluated at the points (x+0.5s.sub.x,y+0.5s.sub.y,t), while the rest of the functions are evaluated at the points (x,y,t) for every (x,y).epsilon.R.sup.2 (h.sub..sigma.,1,h.sub..sigma.,2), t.epsilon.T, s.tbd.(s.sub.x,s.sub.y).epsilon.S. The functions g.sub.t.sup..sigma., g.sub.tu.sup..sigma., g.sub.tv.sup..sigma. are defined, respectively, by the relations ##EQU54## The functions .parallel..gradient.u.sup..sigma. .parallel..sup.2, .parallel..gradient.v.sup..sigma. .parallel..sup.2 are defined as in the relations ##EQU55## The functions (s,.gradient.u.sup..sigma.), (s,.gradient.v.sup..sigma.) are defined by the relations

(s,.gradient.u.sup..sigma.).tbd.(s,.gradient.u.sup..sigma. (x+0.5s.sub.x,y+0.5s.sub.y,t))=.rho..sub.s (u.sup..sigma. (x+s.sub.x,y+s.sub.y,t)-u.sup..sigma. (x,y,t)), (6-34) (s,.gradient.v.sup..sigma.).tbd.(s,.gradient.v.sup..sigma. (x+0.5s.sub.x,y+0.5s.sub.y,t))=.rho..sub.s (v.sup..sigma. (x+s.sub.x,y+s.sub.y,t)-v.sup..sigma. (x,y,t)), (6-35) while the functions (s,.gradient.(.DELTA.u.sup..sigma.)), (s,.gradient.(.DELTA.v.sup..sigma.)) are defined as in the relations (s,.gradient.(.DELTA.u.sup..sigma.)).tbd.(s,.gradient.(.DELTA.u.sup..sigma. (x+0.5s.sub.x,y+0.5s.sub.y,t)))=.rho..sub.s (.DELTA.u.sup..sigma. (x+s.sub.x,y+s.sub.y,t)-.DELTA.u.sup..sigma. (x,y,t)), (6-36) (s,.gradient.(.DELTA.v.sup..sigma.)).tbd.(s,.gradient.(.DELTA.v.sup..sigma. (x+0.5s.sub.x,y+0.5s.sub.y,t)))=.rho..sub.s (.DELTA.v.sup..sigma. (x+s.sub.x,y+s.sub.y,t)-.DELTA.v.sup..sigma. (x,y,t)). (6-37) The function b.sup.2 (s,.gradient.'g.sub.t.sup..sigma.).sup.2 is defined by the relation ##EQU56## where each function (s,.gradient.'g.sub.t.sup..sigma.).sup.2, g.sub.t .epsilon.G.sub.t is determined as (s,.gradient.'g.sub.t.sup..sigma.).sup.2 .tbd.(s,.gradient.'g.sub.t.sup..sigma. (x+0.5s.sub.x,y+0.5s.sub.y,t,u.sup..sigma.,v.sup..sigma.)).sup.2 =(s,.gradient.'g.sub.t.sup..sigma. (x,y,t,u.sup..sigma.,v.sup..sigma.)).sup.2 +(s,.gradient.'g.sub.t.sup..sigma. (x+s.sub.x,y+s.sub.y,t,u.sup..sigma.,v.sup..sigma.)).sup.2. (6-39)

Above, for each value (x',y',t) (equal either to (x,y,t) or to (x+s.sub.x,y+s.sub.y,t)) the value of the functions (s,.gradient.'g.sub.t.sup..sigma. (x',y',t,u.sup..sigma.,v.sup..sigma.)).sup.2 is given by the relation (4-9).

The functions .gamma..sub.i,j.sup..sigma., i, j=1, 2 are defined as in the relation ##EQU57## where .mu. is a small positive constant, while .gamma..sub.i,j.sup..sigma., i, j=1, 2 are defined as in the relations ##EQU58##

The restriction of the system of nonlinear equations (6-27) and the restriction of the system of linear equations (6-28) on the subset .OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2) of the set R.sup.2 (h.sub..sigma.,1,h.sub..sigma.,2) are both accomplished through introduction of the boundary condition as follows. For every image point (x,y).epsilon..OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2) and for every element s.tbd.(s.sub.x,s.sub.y).epsilon.S such that the point (x+s.sub.x,y+s.sub.y) is outside of the subset .OMEGA.(h.sub..sigma.,1,h.sub..sigma., 2) the value of the estimate of the velocity vector (u.sup..sigma. (x+s.sub.x,y+s.sub.y,t),v.sup..sigma. (x+s.sub.x,y+s.sub.y,t) is defined to be identically equal to the value of the estimate of the velocity vector (u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)). Under these conditions the system of nonlinear equations (6-27) and the system of linear equations (6-28) are well defined for every image point (x,y) from the set .OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2) and, therefore, form, respectively, the system of nonlinear equations relative to the unknown estimate of the velocity vector field {(u.sup..sigma. (x,y,t),v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2)} corresponding to the parameter value .sigma. and its finite differences, and the system of linear equations relative to the unknown vector field {(.DELTA.u.sup..sigma. (x,y,t),.DELTA.v.sup..sigma. (x,y,t)).vertline.(x,y).epsilon..OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2)} corresponding to the parameter value .sigma. and its finite differences.

Let us denote by the symbol .vertline.S.vertline. the number of elements of the set S, then for every s.tbd.(s.sub.x,s.sub.y).epsilon.S the constant .rho..sub.s can be defined as in the relation ##EQU59##

The following are two commonly used ways of selecting the set S

S.tbd.{h.sub..sigma.,1,-h.sub..sigma.,1, h.sub..sigma.,2 -h.sub..sigma.,2 }, (6-46) S.tbd.{h.sub..sigma.,1,-h.sub..sigma.,1,h.sub..sigma.,2 -h.sub..sigma.,2,h.sub..sigma.,1 +h.sub..sigma.,2,-h.sub..sigma.,1 -h.sub..sigma.,2,h.sub..sigma.,1 -h.sub..sigma.,2,-h.sub..sigma.,1 +h.sub..sigma.,2 }. (6-47) If the set S is selected as in the relation (6-46) we say that the "five-points finite-difference discretization" is used, while for the selection specified by the relation (6-47) we say that the "nine-points finite-difference discretization" is used.

Let k=0, . . . , n-1 and let {(u.sup..sigma..sbsp.k (x,y,t),v.sup..sigma..sbsp.k (x,y,t)).vertline.(x,y).di-elect cons..OMEGA.(h.sub..sigma..sbsb.k.sub.,1,h.sub..sigma..sbsb.k.sub.,2)} be an estimate of the velocity vector field corresponding to the parameter value .sigma..sub.k. Then the initial estimate {(u.sup..sigma..sbsp.k+1 (x,y,t),v.sup..sigma..sbsp.k+1 (x,y,t)).vertline.(x,y).di-elect cons..OMEGA.(h.sub..sigma..sbsb.k+1.sub.,1,h.sub..sigma..sbsb.k+1.sub.,2)} of the velocity vector field corresponding to the parameter value .sigma..sub.k+1 is defined as follows. For every (x,y).di-elect cons..OMEGA.(h.sub..sigma..sbsb.k+1.sub.,1,h.sub..sigma..sbsb.k+1.sub.,2) such that (x,y).di-elect cons..OMEGA.(h.sub..sigma..sbsb.k.sub.,1,h.sub..sigma..sbsb.k.sub.,2) the value of the estimate (u.sup..sigma..sbsp.k+1 (x,y,t),v.sup..sigma..sbsp.k+1 (x,y,t)) is given by the relation

(u.sup..sigma..sbsp.k+1 (x,y,t),v.sup..sigma..sbsp.k+1 (x,y,t))=(u.sup..sigma..sbsp.k (x,y,t),v.sup..sigma..sbsp.k (x,y,t))(6-48) while for every (x,y).di-elect cons..OMEGA.(h.sub..sigma..sbsb.k+1.sub.,1,h.sub..sigma..sbsb.k+1.sub.,2) such that (x,y).epsilon slash..OMEGA.(h.sub..sigma..sbsb.k.sub.,1,h.sub..sigma..sbsb.k.sub.,2) the value of the estimate (u.sup..sigma..sbsp.k+1 (x,y,t),v.sup..sigma..sbsp.k+1 (x,y,t)) is given as in the relation ##EQU60## where (x.sub.i.sbsb.b.sub.i.sbsb.2.sub., y.sub.i.sbsb.b.sub.i.sbsb.2)=(x,y)+i.sub.1 h.sub..sigma..sbsb.k+1.sub.,1 +i.sub.2 h.sub..sigma..sbsb.k+1.sub.,2. (6-50)

The system of linear equations (6-28) is symmetric and positive definite. As shown below, the system of linear equations (6-28) can be partitioned into two subsystems, one of which dominates the other and can be directly inverted. Based on these properties a fast preconditioned iterative method is constructed, which can be implemented as a fully local, uniform, and parallel algorithm for solving the system of linear equations (6-28). This method is further accelerated in an embodiment of the present invention with the help of the Chebyshev polynomials, and such acceleration retains the above mentioned properties.

If we introduce an ordering on the set .OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2), then the system of linear equations (6-28) can be expressed in the matrix notation as

M.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T =-F.sup..sigma. (u.sup..sigma., v.sup..sigma.), (7-1) where (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T is the transpose vector of the row-vector (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.), corresponding to the vector field {(.DELTA.u.sup..sigma. (x,y,t),.DELTA.v.sup..sigma. (x,y,t)).vertline.(x,y).di-elect cons..OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2)}, which is obtained by first listing the elements of the set {.DELTA.u.sup..sigma. (x,y,t).vertline.(x,y).di-elect cons..OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2)} according to the specified ordering, then listing the elements of the set {.DELTA.v.sup..sigma. (x,y,t).vertline.(x,y).di-elect cons..OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2)} according to the same ordering. The "natural ordering" on the set .OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2) can be defined as follows. The relations (6-2) and (6-2) imply that every element (x,y).di-elect cons..OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2) can be uniquely represented in the form (x,y).tbd.i.sub.1 h.sub..sigma.,1 +i.sub.2 h.sub..sigma.,2 where i.sub.1, and i.sub.2 are integer numbers. Then the element (x,y).tbd.i.sub.1 h.sub..sigma.,1 +i.sub.2 h.sub..sigma.,2 .di-elect cons..OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2) follows the element (x',y').tbd.i'.sub.1 h.sub..sigma.,1 +i'.sub.2 h.sub..sigma.,2 .di-elect cons..OMEGA.(h.sub..sigma.,1,h.sub..sigma.,2) if i.sub.2 >i'.sub.2 or if i.sub.2 =i'.sub.2 and i.sub.1 >i'.sub.1.

By partitioning the square matrix M.sup..sigma. (u.sup..sigma., v.sup..sigma.) into four square submatrices corresponding to the vectors (.DELTA.u.sup..sigma.) and (.DELTA.v.sup..sigma.), respectively, and by partitioning the vector -F.sup..sigma. (u.sup..sigma., v.sup..sigma.) into two subvectors the relation (7-1) representing the system of linear equations (6-28) becomes ##EQU61##

Let the matrix D.sub.i,j.sup..sigma. (u.sup..sigma.,v.sup..sigma.) be the diagonal part of the matrix M.sub.i,j.sup..sigma. (u.sup..sigma.,v.sup..sigma.), while the matrix -B.sub.i,j.sup..sigma. (u.sup..sigma., v.sup..sigma.) is the off-diagonal part of the matrix M.sub.i,j.sup..sigma. (u.sup..sigma., v.sup..sigma.), so that the relation

M.sub.i,j.sup..sigma. (u.sup..sigma.,v.sup..sigma.)=D.sub.i,j.sup..sigma. (u.sup..sigma.,v.sup..sigma.)-B.sub.i,j.sup..sigma. (u.sup..sigma.,v.sup..sigma.) (7-4) is satisfied for every i, j=1, 2, and let the matrices D.sup..sigma. (u.sup..sigma.,v.sup..sigma.), B.sup..sigma. (u.sup..sigma.,v.sup..sigma.) be defined as in the relations ##EQU62## then the system of linear equations (6-28) takes the form (D.sup..sigma. (u.sup..sigma.,v.sup..sigma.)-B.sup..sigma. (u.sup..sigma.,v.sup..sigma.))(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T =-F.sup..sigma. (u.sup..sigma.,v.sup..sigma.). (7-6)

The quadratic form corresponding to the system of linear equations (6-28) represented by the relation (7-6) can be given as in the relation

(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)(D.sup..sigma. (u.sup..sigma., v.sup..sigma.)-B.sup..sigma. (u.sup..sigma.,v.sup..sigma.))(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T (7-7) and is greater than the quadratic form ##EQU63## If the relations (s,.gradient.(.DELTA.u.sup..sigma.)).tbd.(s,.gradient.(.DELTA.u.sup..sigma. (x+0.5s.sub.x,y+0.5s.sub.y,t)))=-.rho..sub.s (.DELTA.u.sup..sigma. (x+s.sub.x,y+s.sub.y,t)+.DELTA.u.sup..sigma. (x,y,t)), (7-9) (s,.gradient.(.DELTA.v.sup..sigma.)).tbd.(s,.gradient.(.DELTA.v.sup..sigma. (x+0.5s.sub.x,y+0.5s.sub.y,t)))=-.rho..sub.s (.DELTA.v.sup..sigma. (x+s.sub.x,y+s.sub.y,t)+.DELTA.v.sup..sigma. (x,y,t)), (7-10) are used instead of the relations (6-36), (6-37) to define the functions (s,.gradient.(.DELTA.u.sup..sigma.)), (s,.gradient.(.DELTA.v.sup..sigma.)) appearing in the system of linear equations (6-28); then the system of linear equations (6-28) takes the form D.sup..sigma. (u.sup..sigma.,v.sup..sigma.)+B.sup..sigma. (u.sup..sigma.,v.sup..sigma.))(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T =-F.sup..sigma. (u.sup..sigma.,v.sup..sigma.). (7-11) The quadratic form corresponding to the system of linear equations (6-28) represented by the relation (7-11) can be given as in the relation (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)(D.sup..sigma. (u.sup..sigma.,v.sup..sigma.)+B.sup..sigma. (u.sup..sigma.,v.sup..sigma.))(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T (7-12) and is greater than the quadratic form (7-8).

It is not difficult to see that the quadratic form (7-8) is equal to the quadratic form

(1-.lambda.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)D.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T, (7-13) where the positive constant .lambda. is given by the relation .lambda.=1/(.mu.+1). (7-14)

By adding up the expressions of the relations (7-7) and (7-12) it can be seen that the matrix D.sup..sigma. (u.sup..sigma.,v.sup..sigma.) is positive definite which implies that the diagonal matrix

D.sub.1,1.sup..sigma. (u.sup..sigma.,v.sup..sigma.)D.sub.2,2.sup..sigma. (u.sup..sigma.,v.sup..sigma.)-D.sub.1,2.sup..sigma. (u.sup..sigma.,v.sup..sigma.)D.sub.2,1.sup..sigma. (u.sup..sigma.,v.sup..sigma.) (7-15) as well as diagonal matrices D.sub.1,1.sup..sigma. (u.sup..sigma.,v.sup..sigma.), D.sub.2,2.sup..sigma. (u.sup..sigma.,v.sup..sigma.) are positive definite.

Let C.sup..sigma. (u.sup..sigma.,v.sup..sigma.) be a block lower triangular matrix defined by the relations ##EQU64## then the following relation is satisfied

C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.T C.sup..sigma. (u.sup..sigma.,v.sup..sigma.)=D.sup..sigma. (u.sup..sigma.,v.sup..sigma.).(7-18)

Let (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) be a vector defined by the relation

(.DELTA.u.sup..sigma.,v.sup..sigma.).sup.T =C.sup..sigma.)u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T. (7-19)

If both sides of the system of linear equations (7-6) of the unknown vector (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) are multiplied on the matrix C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-T from the left and if the unknown vector (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) is substituted with the unknown vector (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) the system of linear equations (7-6) becomes the following system of linear equations relative to the unknown vector (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.):

(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T =C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-T B.sup..sigma. (u.sup..sigma.,v.sup..sigma.)C.sup..sigma. (i u.sup..sigma.,v.sup..sigma.).sup.-1 (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T -C.sup..sigma. (u.sup..sigma., v.sup..sigma.).sup.-T F.sup..sigma. (u.sup..sigma.,v.sup..sigma.). (7-20)

By substitution of the unknown vector (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) with the unknown vector (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) the quadratic form (7-7) becomes the quadratic form

(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T -(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-T B.sup..sigma. (u.sup..sigma.,v.sup..sigma.)C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-1 (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T, (7-21) the quadratic form (7-12) becomes the quadratic form (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T +(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-T B.sup..sigma. (u.sup..sigma.,v.sup..sigma.)C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-1 (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T, (7-22) and the quadratic form (7-13) becomes (1-.lambda.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T. (7-23)

The fact that the quadratic forms (7-21) and (7-22) are both greater than the quadratic form (7-23) implies that the following relations are satisfied:

(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-T B.sup..sigma. (u.sup..sigma.,v.sup..sigma.)C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-1 (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T <+.lambda.(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T, (7-24) (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-T B.sup..sigma. (u.sup..sigma.,v.sup..sigma.)C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-1 (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T >-.lambda.(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.)(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T. (7-25)

The relations (7-24) and (7-25), in turn, imply that the L.sup.2 norm of the matrix C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-T B.sup..sigma. (u.sup..sigma.,v.sup..sigma.)C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-1 is less than .lambda., i.e.,

.parallel.C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-T B.sup..sigma. (u.sup..sigma.,v.sup..sigma.)C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-1 .parallel..sub.2 <.lambda.<1.(7-26)

Let G.sup..sigma. (u.sup..sigma.,v.sup..sigma.) be a matrix defined as in the relation

G.sup..sigma. (u.sup..sigma.,v.sup..sigma.)=C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-T B.sup..sigma. (u.sup..sigma.,v.sup..sigma.) C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-1, (7-27) while h.sup..sigma. (u.sup..sigma.,v.sup..sigma.) is a vector given by the relation h.sup..sigma. (u.sup..sigma.,v.sup..sigma.)=C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-T F.sup..sigma. (u.sup..sigma.,v.sup..sigma.); (7-28) then the following basic iterative method can be used to obtain an approximation to the solution of the system of linear equations (6-28).

The initial approximation (.DELTA.u.sub.0.sup..sigma.,.DELTA.v.sub.0.sup..sigma.) to the solution of the system of linear equations (7-20) is defined to be identically equal to zero. For every n=0,1, . . . the approximation (.DELTA.u.sub.n+1.sup..sigma.,.DELTA.v.sub.n+1.sup..sigma.) is defined in terms of the approximation (.DELTA.u.sub.n.sup..sigma.,.DELTA.v.sub.n.sup..sigma.) as in the relation

(.DELTA.u.sub.n+1.sup..sigma.,.DELTA.v.sub.n+1.sup..sigma.).sup.T =G.sup..sigma. (u.sup..sigma.,v.sup..sigma.)(.DELTA.u.sub.n.sup..sigma.,.DELTA.v.sub.n.sup..sigma.).sup.T h.sup..sigma. (u.sup..sigma.,v.sup..sigma.).(7-29)

The process is continued until a proper approximation (.DELTA.u.sub.N.sup..sigma.,.DELTA.v.sub.N.sup..sigma.) to the solution of the system of linear equations (7-20) is achieved. The approximation (.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.) to the solution of the system of linear equations (6-28) is defined by the relation

(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T =C.sup..sigma. (u.sup..sigma.,v.sup..sigma.).sup.-1 (.DELTA.u.sub.N.sup..sigma.,v.sub.N.sup..sigma.).sup.T. (7-30)

The performance of the basic iterative method (7-29) can be improved with the help of the polynomial acceleration applied to the basic iterative method. Two methods are described below: a conjugate gradient polynomial acceleration, and a Chebyshev polynomial acceleration. For the sake of simplicity of the presentation the following notation shall be used:

w.tbd.(.DELTA.u.sup..sigma.,.DELTA.v.sup..sigma.).sup.T, G.tbd.G.sup..sigma. (u.sup..sigma.,v.sup..sigma.), h.tbd.h.sup..sigma. (u.sup..sigma.,v.sup..sigma.). (7-31) Then the basic iterative method (7-29) becomes w.sub.n+1 =Gw.sub.n -h. (7-32)

The polynomial acceleration of the basic iterative method (7-32) can be described as follows. Start with the initial approximation w.sub.0 which is identically equal to zero, then for every n=0,1, . . . the following iterative procedure is applied:

w.sub.n+1 =.rho..sub.n (.gamma..sub.n (Gw.sub.n -h)+(1-.gamma..sub.n)w.sub.n)+(1-.rho..sub.n)w.sub.n-1. (7-33) In the case of the conjugate gradient method the coefficients .rho..sub.n and .gamma..sub.n are given by the relation ##EQU65## while in the case of the Chebyshev method the coefficients .rho..sub.n, and .gamma..sub.n are given by the relation ##EQU66##

The conjugate gradient method can be expressed in the alternative form

w.sub.n+1 +w.sub.n +.alpha..sub.n p.sub.n, r.sub.n+1 =r.sub.n -.alpha..sub.n q.sub.n, (7-39) where the coefficients .alpha..sub.n are defined as in the relation (7-35), while the vectors p.sub.n, q.sub.n are given by the relation (7-36). The alternative form for the Chebyshev method is defined as w.sub.n+1 =w.sub.n +p.sub.n, r.sub.n+1 =r.sub.n -q.sub.n, (7-40) where p.sub.n =.rho..sub.n r.sub.n +(.rho..sub.n -1)p.sub.n-1, n.gtoreq.1, p.sub.0 =r.sub.0, q.sub.n =p.sub.n -Gp.sub.n, n.gtoreq.0, (7-41) and the coefficients .rho..sub.n are given by the relation (7-38).

In an exemplary embodiment of the invention, the five level resolution pyramid was used with the value of the parameter .sigma. decreasing by a factor of 2.0 for each successively finer resolution level. A nine-points finite-difference discretization was used. The positive integer constant .kappa.' appearing in the relations (6-6), (6-5) is equal to 4, while the positive integer constant .kappa." appearing in the relations (6-19), (6-20) is equal to 2. With these values, the present invention is capable of estimating the velocity vector field with good accuracy in the case of the images of real world scenes containing multiple occluding boundaries and objects undergoing nonrigid body motions.

The mathematical analysis and methods discussed above are used to modify a time-varying image sequence in the embodiments of the method of the present invention described below and illustrated in FIGS. 2-8.

FIG. 2 illustrates an embodiment of the method of the present invention for modifying at least one image in a time-varying image sequence. In step 40, a time-varying image sequence is input. In step 42, optical flow constraints are combined with directional smoothness constraints and regularization constraints in a functional of the estimate of the parametric velocity vector field such that weighting coefficients of the functional are functions of the parametric velocity vector field to be estimated. Then, in step 44, a finite sequence of parameter values are selected for solving a system of nonlinear equations that arise from optimality criterion of the functional, wherein the solution of the system of nonlinear equations is smooth with respect to a parameter value. Then, in step 46, the system of nonlinear equations is successively solved starting with a parameter value where the solution of the system of nonlinear equations is equal to a given vector constant and proceeds in a direction decreasing with respect to the parameter values using the solution of the system of nonlinear equations obtained for the preceding parameter value as the initial estimate. As shown in step 48, the estimate of the parametric velocity vector field is applied in order to modify at least one image in the image sequence.

FIG. 3 shows a different embodiment of the method of the present invention in which the first two steps, steps 50 and 52 are the same as steps 40 and 42 in the embodiment of FIG. 2. In step 54, a finite sequence of parameter values is selected for solving a system of nonlinear equations that arise from optimality criterion applied to the functional. The system of nonlinear equations is such that a quadratic form associated with a Jacobian of no operator forming the system of nonlinear equations is not positive definite near the occluding boundaries and is asymmetric.

An initial estimate of the velocity vector field is selected corresponding to the initial parameter value and takes the initial parameter value as the current parameter value in step 56. Then, in step 58, the system of nonlinear equations is formed corresponding to the current parameter value and the initial estimate of the velocity vector field is taken as the current estimate of the velocity vector field. In step 60, the Jacobian of the operator forming a system of nonlinear equations is approximated with a symmetric and positive definite linear operator. This linear operator is used to form a system of linear equations.

In step 62, a preconditioned iterative process is performed for solving the system of linear equations by either a conjugate gradient acceleration or a Chebyshev polynomial acceleration, and then the solution is used to improve the current estimate. Step 64 is a decision step, in which it is determined whether the system of nonlinear equations is solved. If the determination is negative, the flow loops back to between steps 58 and 60. If the determination is positive, the method proceeds to step 66, a second decision step, in which it is determined whether the current parameter value is final. If this determination is no, the current parameter value is incremented to then next parameter value in step 68. The estimate of the velocity vector field corresponding to the current parameter value is projected to the estimate of the velocity vector field corresponding to the next parameter value. The flow then loops back to between steps 56 and 58.

If the current parameter value is final, as determined in step 66, then the estimate of the parameter velocity vector field is applied in step 70 to modify at least one image of the time-varying image sequence.

FIG. 4 illustrates another embodiment of the present invention in which the first two steps, steps 72 and 74, are the same as steps 40 and 42 of the embodiment of FIG. 3. In step 76, a system of nonlinear equations is solved which arise from optimality criterion of the functional by selecting a finite sequence of parameter values and successively solving the system of linear equations. These linear equations are successively solved starting with a parameter value where the solution of the system of nonlinear equations is equal to a given vector constant and proceed in a direction decreasing with respect to the parameter values using the solution of the system of nonlinear equations obtained for the preceding parameter value as the initial estimate. This includes performing a finite-difference discretization of the system of nonlinear equations. Then, in step 78, the estimate of the parametric velocity vector field is applied to modify at least one image in the image sequence.

In FIG. 5, another embodiment of the method of the present invention for modifying an image in image sequence is illustrated. In step 80, a time-varying image sequence is input. Optical flow constraints, directional smoothness constraints and regularization constraints are combined in step 82 in a functional of the estimate of the parametric velocity vector field such that weighting coefficients of the functional are functions of the parametric velocity vector field to be estimated. A dual feed-back interaction is provided between the optical flow constraints and the directional smoothness constraints including selectively applying variation to the optical flow constraints and to the directional smoothness constraints and not to the weighting coefficients. This permits using unknown estimates of the velocity vectors as part of variables in the weighting coefficients.

In step 84, a system of nonlinear equations is formed that arise from optimality criterion of the functional, where the solution of nonlinear equations is smooth with respect to a parameter value. The system of nonlinear equations is successively solved in step 86, starting with a parameter value where the solution of the system of nonlinear equations is equal to a given vector constant and proceeds in a direction decreasing with respect to the parameter values using the solution of the system of nonlinear equations obtained for the preceding parameter value as the initial estimate. Then, in step 88, the estimate of the parametric velocity vector field is applied to modify at least one image in the image sequence.

FIG. 6 illustrates a flow chart of another embodiment of the method of the present invention in which a time-varying image sequence is input in step 90. In step 92, optical flow constraints, directional smoothness constraints and regularization constraints are combined in a functional of the estimate of the parametric velocity vector field such that weighting coefficients of the functional are functions of the parametric velocity vector field to be estimated, and a dual feed-back interaction is provided between the optical flow constraints and the directional smoothness constraints.

In step 94, a finite sequence of parameter values is selected for solving a system of nonlinear equations that arise from optimality criterion applied to the functional. This system of nonlinear equations is such that a quadratic form associated with a Jacobian of an operator forming the system of nonlinear equations is not positive definite near the occluding boundaries and is asymmetric. In step 96, an initial estimate of the velocity vector field corresponding to the initial parameter value is selected, and the initial parameter value is taken as the current parameter value.

The system of nonlinear equations corresponding to the current parameter value is formed in step 98, and the initial estimate of the velocity vector field is the current estimate of the velocity vector field as taken.

In step 100, the Jacobian of the operator forming the system of nonlinear equations is approximated with a symmetric and positive definite linear operator. This linear operator is used to form a system of linear equations. In step 102, a preconditioned iterative process is performed for solving the system of linear equations by either a conjugate gradient acceleration or a Chebyshev polynomial acceleration. This solution is then used to improve the current estimate.

It is determined in step 104 whether the system of nonlinear equations is solved. If the determination is no, the flow loops back to between steps 98 and 100. If the determination is yes, a second decision step, step 106, is reached in which it is determined whether the current parameter value is final. If this determination is no, the current parameter value is incremented in step 108 to the next parameter value. The estimate of the velocity vector field to the current parameter value is projected to the estimate of the velocity vector field corresponding to the next parameter value. The flow then loops back to between steps 96 and 98.

If the current parameter value is final, as determined in step 106, the method proceeds to step 110, in which the estimate of the parametric velocity vector field is applied in order to modify at least one image of the time-varying image sequence.

In the embodiment of FIG. 7, steps 112, and 114 are the same as steps 90, and 92 in FIG. 6. In step 116, a system of nonlinear equations is formed that arise from optimality criterion of the functional. The solution of nonlinear equations is smooth with respect to a parameter value.

In step 118, the system of nonlinear equations is successively solved starting with a parameter value where the solution of the system of nonlinear equations is equal to a given vector constant and proceeds in the direction decreasing with respect to the parameter values using the solution of the system of nonlinear equations obtained for the preceding parameter value as the initial estimate. This includes performing a finite-difference discretization of the system of nonlinear equations. Then, in step 120, the estimate of the parametric velocity vector field is applied to modify at least one image in the image sequence.

FIG. 8 illustrates another embodiment of the method of the present invention. In step 122, the images of a time-varying image sequence are digitized. Successive digitized images of the image sequence are treated as generalized functions in step 124. A specific subset of a set of testing functions is obtained in step 126 by introducing a scaling factor into one of the testing functions that is a measurement function. The derivatives of the digitized images are treated in step 128 as generalized derivatives observed on the specific subset of the set of testing functions. In step 130, a system of nonlinear equations corresponding to the testing functions are sequentially solved with a progressively lower degree of measurement blur. Finally, in step 132, the estimate of the parametric velocity vector field is applied to modifying at least one image in the image sequence.

While there has been shown what is considered to be the preferred embodiments of the invention, it will be manifest that many changes and modifications may be made therein without departing from the essential spirit of the invention. It is intended, therefore, in the annexed claims, to cover all such changes and modifications as may fall within the true scope of the invention.

Claims

1. A method of modifying a time-varying image sequence comprising the steps of:

a. estimating a parametric velocity vector field that characterizes changes in successive images of the image sequence, said step of estimating including:

i. combining optical flow constraints, directional smoothness constraints and regularization constraints in a functional of the estimate of the parametric velocity vector field such that weighting coefficients of the functional are functions of the parametric velocity vector field to be estimated; and

ii. solving a system of nonlinear equations that arise from optimality criterion of the functional;

b. applying the estimate of the parametric velocity vector field to modify at least one image in the image sequence.

2. The method of claim 1, wherein the step of solving the system of nonlinear equations includes the steps of selecting a finite sequence of parameter values, and successively solving the system of nonlinear equations starting with a parameter value where the solution of the system of nonlinear equations is equal to a given vector constant and proceeding in a direction decreasing with respect to the parameter values using the solution of the system of nonlinear equations obtained for the preceding parameter value as the initial estimate.

3. The method of claim 1, wherein the solution of the system of nonlinear equations is smooth with respect to a parameter value.

4. The method of claim 1, wherein the system of nonlinear equations is such that a quadratic form associated with a Jacobian of an operator forming the system of nonlinear equations is not positive definite near the occluding boundaries and is asymmetric.

5. The method of claim 4, further comprising approximating the Jacobian with a symmetric and positive definite linear operator.

6. The method of claim 5, further comprising deriving a parametric system of linear equations from the Jacobian.

7. The method of claim 6, further comprising performing a preconditioned iterative process for solving the system of linear equations.

8. The method of claim 7, wherein the iterative process is a conjugate gradient acceleration.

9. The method of claim 7, wherein the iterative process is a Chebyshev polynomial acceleration.

10. The method of claim 2, wherein the successive solving of the system of nonlinear equations includes performing a finite-difference discretization of the system of nonlinear equations.

11. The method of claim 1, wherein the system of nonlinear equations that are solved is: ##EQU67## where the parameters appearing in this system of nonlinear equations variables are defined as follows: .sigma..delta.(.sigma.) are smoothing parameters; u.sup..sigma. is a horizontal component of the estimate of the velocity vector, v.sup..sigma. is a vertical component of the estimate of the velocity vector; G.sub.1 is a set of optical flow constraints; r, p, q, .rho.(g.sub.1), g.sub.1 .di-elect cons.G.sub.1 are constant parameters specifying the part of the system of nonlinear equations representing the optical flow constraints; .parallel..gradient.u.sup..sigma. .parallel. is a norm of the gradient of the horizontal component of the estimate of the velocity vector; .parallel..gradient.v.sup..sigma. .parallel. is a norm of the gradient of the vertical component of the estimate of the velocity vector, g.sub.1.sup..sigma., g.sub.1 .di-elect cons.G.sub.1 is an optical flow constraint; g.sub.1u.sup..sigma., g.sub.1 .di-elect cons.G.sub.1 is a derivative of the optical flow constraint g.sub.1.sup..sigma. with respect to the horizontal component u.sup..sigma. of the estimate of the velocity vector; g.sub.1v.sup..sigma., g.sub.1 .di-elect cons.G.sub.1 is a derivative of the optical flow constraint g.sub.1.sup..sigma. with respect to the vertical component v.sup..sigma. of the estimate of the velocity vector; S is a set of directions used in the directional smoothness constraints; a,c,b,.rho..sub.s,s.di-elect cons.S are constant parameters specifying the part of the system of nonlinear equations representing the directional smoothness constraints; b.sup.2 (s,.sigma.'g.sub.1.sup..sigma.).sup.2 is a weighted average of the squares of the directional derivatives of the optical flow constraints in the direction s; (s.gradient.u.sup..sigma.) is a directional derivative of the horizontal component u.sup..sigma. of the estimate of the velocity vector in the direction s; (s,.gradient.v.sup..sigma.) is a directional derivative of the vertical component v.sup..sigma. of the estimate of the velocity vector in the direction s; .gamma..sup..sigma. is the constant parameter.

12. A method of modifying a time-varying image sequence comprising the steps of:

a. estimating a parametric velocity vector field that characterizes changes in successive images of the image sequence, said step of estimating including:

i. combining optical flow constraints, directional smoothness constraints and regularization constraints in a functional of the estimate of the parametric velocity vector field such that weighting coefficients of the functional are functions of the parametric velocity vector field to be estimated;

ii. providing a dual feed-back interaction between the optical flow constraints and the directional smoothness constraints;

iii. solving a system of nonlinear equations that arise from optimality criterion of the functional;

b. applying the estimate of the parametric velocity vector field to modify at least one image in the image sequence.

13. The method of claim 12, wherein the step of providing dual feed-back interaction includes selectively applying variation to the optical flow constraints and to the directional smoothness constraints and not to the weighting coefficients, to thereby permit using unknown estimates of the velocity vectors as part of variable in the weighting coefficients.

14. The method of claim 12, wherein the step of solving the system of nonlinear equations includes the steps of selecting a finite sequence of parameter values, and successively solving the system of nonlinear equations starting with a parameter value where the solution of the system of nonlinear equations is equal to a given vector constant and proceeding in a direction decreasing with respect to the parameter values using the solution of the system of nonlinear equations obtained for the preceding parameter value as the initial estimate.

15. The method of claim 12, wherein the solution of the system of nonlinear equations is smooth with respect to a parameter value.

16. The method of claim 12, wherein the system of nonlinear equations is such that a quadratic form associated with a Jacobian of an operator forming the system of nonlinear equations is not positive definite near the occluding boundaries and is asymmetric.

17. The method of claim 16, further comprising approximating the Jacobian with a symmetric and positive definite linear operator.

18. The method of claim 17, further comprising deriving a parametric system of linear equations from the Jacobian.

19. The method of claim 18, further comprising performing a preconditioned iterative process for solving the system of linear equations.

20. The method of claim 19, wherein the iterative process is a conjugate gradient acceleration.

21. The method of claim 19, wherein the iterative process is a Chebyshev polynomial acceleration.

22. The method of claim 14, wherein the successive solving of the system of nonlinear equations includes performing a finite-difference discretization of the system of nonlinear equations.

23. The method of claim 12, wherein the system of nonlinear equations that are solved is: ##EQU68## where the parameters appearing in this system of nonlinear equations variables are defined as follows: .sigma..delta.(.sigma.) are smoothing parameters; u.sup..sigma. is a horizontal component of the estimate of the velocity vector, v.sup..sigma. is a vertical component of the estimate of the velocity vector; G.sub.1 is a set of optical flow constraints; r, p, q, .rho.(g.sub.1), g.sub.1 .di-elect cons.G.sub.1 are constant parameters specifying the part of the system of nonlinear equations representing the optical flow constraints; .parallel..gradient.u.sup..sigma. .parallel. is a norm of the gradient of the horizontal component of the estimate of the velocity vector; .parallel..gradient.v.sup..sigma. .parallel. is a norm of the gradient of the vertical component of the estimate of the velocity vector, g.sub.1.sup..sigma., g.sub.1 .di-elect cons.G.sub.1 is an optical flow constraint; g.sub.1u.sup..sigma., g.sub.1 .di-elect cons.G.sub.1 is a derivative of the optical flow constraint g.sub.1.sup..sigma. with respect to the horizontal component u.sup..sigma. of the estimate of the velocity vector; g.sub.1v.sup..sigma., g.sub.1 .di-elect cons.G.sub.1 is a derivative of the optical flow constraint g.sub.1.sup..sigma. with respect to the vertical component v.sup..sigma. of the estimate of the velocity vector; S is a set of directions used in the directional smoothness constraints; a,c,b,.rho..sub.s,s.di-elect cons.S are constant parameters specifying the part of the system of nonlinear equations representing the directional smoothness constraints; b.sup.2 (s,.sigma.'g.sub.1.sup..sigma.).sup.2 is a weighted average of the squares of the directional derivatives of the optical flow constraints in the direction s; (s.gradient.u.sup..sigma.) is a directional derivative of the horizontal component u.sup..sigma. of the estimate of the velocity vector in the direction s; (s,.gradient.v.sup..sigma.) is a directional derivative of the vertical component v.sup..sigma. of the estimate of the velocity vector in the direction s; .gamma..sup..sigma. is the constant parameter.

24. A method of modifying a time-varying image sequence comprising the steps of:

a. digitizing images of the time-varying image sequence;

b. estimating a parametric velocity vector field that characterizes changes in successive digitized images of the image sequence, said step of estimating including:

i. treating the digitized images as generalized functions;

ii. obtaining a specific subset of a set of testing functions by introduction of a scaling factor into one of the testing functions that is a measurement function;

iii. treating derivatives of the digitized images as generalized derivatives observed on the specific subset of the set of testing functions; and

iv. sequentially solving with a progressively lower degree of measurement blue a system of nonlinear equations corresponding to the testing functions;

b. applying the estimate of the parametric velocity vector field to modify at least one image in the image sequence.