Segmentation-Based Method for Motion-Compensated Frame Interpolation

The invention relates generally to the field of digital image processing, and in particular to motion compensation in digital motion image processing. More specifically, the invention relates to a method for motion compensation useful in video frame interpolation.

Motion compensation is widely used in digital video processing algorithms to exploit temporal redundancy in an effort to produce high-quality results. An example of video processing is frame interpolation, used for standards conversion as well as for generation of slow and fast motion in motion picture film/video post-production. In frame interpolation, direct temporal interpolation, i.e. without using motion information, results either in blurring or temporal aliasing when there is relative motion between objects from frame-to-frame. Motion-compensated interpolation, on the other hand, is capable of producing results without blur and aliasing.

The success of motion-compensated interpolation hinges on the accuracy of the motion information. Because motion estimation is an underdetermined problem, it relies on a priori constraints and information about the motion vector field (e.g. smoothness), the image intensity, and models for the 3-D-to-2-D projection. When such constraints becomes invalid, motion estimation becomes inaccurate. When a background is occluded with an object moving across the background, fundamental constraints such as image intensity conservation (the optical flow constraint), used in motion estimation, becomes invalid. At the boundaries of the occluded regions, corresponding to object boundaries, smoothness of the motion field is invalid since motion of occluding objects may be significantly different. Occlusions manifest themselves as covered and uncovered regions, or self occlusions due to local deformation in non-rigid objects. Of particular interest for the subject invention are occlusions that are due to covered/uncovered regions, assuming that object deformations are sufficiently small. A depiction of such regions for a hypothetical 1-D case is given in FIGS. 1A and 1B.

FIG. 1A is a depiction of an occlusion in a 1-D image. As the rectangular object 10 in Frame 1 at time t.sub.P, moves against a stationary background 12 to a new location indicated by the location of the rectangular object 10' in Frame 2, at time t.sub.Q, FIG. 1B, it covers region 16 and uncovers region 14. The uncovered region 14 is absent at time t.sub.P and the covered region 16 is absent at time t.sub.Q. The 1-D example is shown for simplicity with the understanding that 2-D and 3-D occlusions may also occur.

FIG. 2 illustrates an Intermediate Frame, occurring at time t, wherein the region 16 is partially covered and the region 14 is partially uncovered by the rectangular object 10". In other words, the object 10 moved from its position in Frame 1 to the position shown in the Intermediate Frame to its final position shown in Frame 2. Applying linear trajectories 20 between like points at times t.sub.P and t.sub.Q. provides information about the image at time t. To be noted: If the time period between t.sub.P and t.sub.Q is great, any non-linear motion of the rectangular object 10 at time t will not be accurately detected resulting in an inaccurate Intermediate Frame determination. In a like manner, if the object 10 linearly moved to, for example, the far left and linearly back to the position shown in Frame 2, the Intermediate Frame would not represent that motion. Therefore, the time between t.sub.P and t.sub.Q must be relatively small for accurate results when using the present algorithm. Applying motion-compensated, bi-directional interpolation along the linear motion trajectories 20 results in occluded regions 18 in the intermediate frame, denoted by ellipses. Image intensities for the occluded regions are copied from the appropriate existing frame. That is, for the occluded region 18, the image intensities are copied from Frame 2, and for the occluded region 18' image intensities are copied from Frame 1. The trajectories 20 shown are ideal trajectories assuming exact motion estimation.

In summary then, using motion-compensated frame interpolation, a frame at an intermediate time instance, t, is interpolated from existing frames that are captured at times t.sub.P and t.sub.Q, where t.sub.P <t<t.sub.Q. The linear motion trajectories 20 (vectors) are estimated between the existing Frames 1 and 2 such that they describe linear trajectories passing through pixel locations at the missing Intermediate Frame. FIG. 2 assumes the ideal case where motion vectors are accurate and the occluded regions 18 are accurately detected and handled. For linear motion trajectories (vectors) corresponding to non-occluded regions a linear combination (e.g. a weighted average) of the two image values at times t.sub.P and t.sub.Q is used to estimate the missing pixel value through which the associated vector passes. In uncovered regions of the Intermediate Frame, background intensity values from the future frame, Frame 2, are taken as the estimates for the missing pixels. Conversely for covered regions the background intensity values for the past frame Frame 1 are taken as estimates for the missing pixel values. This is what is meant by "occlusion handling" in the context of frame interpolation.

The accuracy of motion estimation and detection of occlusion regions are important practical problems that have been addressed. See, for example, Philippe Robert, "motion compensating interpolation considering occluded, appearing and disappearing areas", Signal Processing of HDTV, III, H. Yasuda and Chiariglione (Eds.), pp. 329-341, Elsevier Publishers, Amsterdam, 1991; and M. Bierling and R. Thoma, "Motion-compensating interpolation considering covered and uncovered background", Signal Processing: Image Communication, pp. 191-212, 1989. In addition, inaccuracies in motion estimation are not limited to only occluded regions. In almost all motion estimation algorithms, the problems within the occluded region propagate to regions outside the occluded region, adversely affecting the estimates immediately outside the occluded region. This case is depicted in FIG. 3. The case of FIG. 3 should be viewed as a non-ideal version of the case shown in FIG. 2 and represents what happens in practice. The regions 22 and 22' containing the occluded regions 18 and 18', respectively, at the Intermediate Frame location of the rectangular object 10" over which motion estimates are inaccurate are designated "bad regions," respectively. In practice, motion estimates are inaccurate within "bad regions", containing the occluded regions 18 and 18' and the immediate neighborhood.

The origin of "bad regions" can be explained as follows: At occluded region boundaries, the motion vector field should not be smooth, however the optical flow constraint is valid. On the other hand, within the occluded region, the optical flow constraint is invalid, but the motion vector field is smooth. Therefore, a successful algorithm has to enforce/suppress these two constraints at the appropriate locations at and around the occluded regions. Most algorithms are incapable of doing so since they tend to enforce these constraints with a fixed relative strength for the entire image region. The Fogel algorithm (see S. Fogel, "The estimation of velocity vector fields from time-varying image sequences, Computer Vision Graphic and Image Processing: Image Understanding, vol. 53, pp. 253-287, 1991), on the other hand, deals with occlusions using an adaptive mechanism that attempts to detect occluded regions and to adjust the relative strengths of these two constraints. Detection of occluded regions is based on implicit detection of image intensity transitions (i.e. edges) and gradients of the motion vector field. Because detection of occluded regions is itself a difficult problem, the Fogel algorithm may be unable, for instance, to turn off the smoothness.sub.-- constraint at the exact position of the occlusion boundary. In fact, in certain cases, detection of edges may be rather difficult for the eye of an observer. This causes the propagation of estimation errors into what has been called the "bad regions" of FIG. 3.

There is a need therefore for an object-based method to accurately determine occluded regions, to remove the bad regions, and to perform appropriate interpolation of the Intermediate Field at occluded regions.

The present invention is directed to overcoming one or more of the problems set forth above. This is achieved by identifying and processing occluding objects separately. The contours of occluding objects (occlusion boundaries) are manually determined via an interactive drawing tool that has a graphical user interface. This segmentation in effect amounts to manually providing the Fogel algorithm with the information about where to enforce/suppress optical flow and smoothness constraints. Briefly summarized, according to one aspect of the present invention, a method of motion compensated frame interpolation wherein an intermediate image frame is generated between two image frames in a motion image sequence, includes the steps of: identifying a foreground and a background for two successive frames in a motion image sequence; estimating a velocity vector field for the foreground, and a velocity vector field for the background; using the estimated velocity vector field for the foreground to identify the location of the foreground and the background in the intermediate image; and obtaining the intermediate image by using the foreground velocity vector field to perform motion compensated interpolation for the foreground in the intermediate image and using the background velocity vector field to perform motion compensated interpolation for the background in the intermediate image.

These and other aspects, objects, features and advantages of the present invention will be more clearly understood and appreciated from a review of the following detailed description of the preferred embodiments and appended claims, and by reference to the accompanying drawings.

The present invention is more robust in the presence of noise in the original image and renders a more accurate intermediate image than Prior Art methods.

FIGS. 1A and 1B are diagrams useful in describing the present invention, showing two successive image frames having a moving object;

FIG. 2 is a diagram similar to that shown in FIGS. 1A and 1B, illustrating the generation of an Intermediate Frame at a time t and the resulting occluded regions;

FIG. 3 is a diagram similar to that shown in FIG. 2, illustrating the problem of estimation of velocity vectors in the occluded regions;

FIG. 4 is a block diagram illustrating an image processing system useful for practicing the present invention;

FIG. 5A and 5B form a flow chart that is useful in understanding the method of the present invention; and

FIGS. 6A and 6B illustrate in diagram form a truck image moving over a background image and the results of that motion is a 1D image plane.

To facilitate understanding, identical reference numerals have been used, where possible, to designate identical elements that are common to the figures.

The method of the present invention is useful as a tool to solve the problems that arise in "difficult" motion clips where errors due to inaccurate detection and mishandling of occlusion are most objectionable (e.g. close ups of faces, etc.). In many other cases, such as detailed scenes viewed from a distance (e.g. an explosion scene) one may not need to resort to this method. It should also be noted that this approach is not intended to solve the problems due to severe self occlusions (e.g. due to local deformations) that may take place within objects rather than at their boundaries.

In what follows, the main principle of the method of the present invention is described. Next, a detailed mathematical description of the method and its implementation is described.

The main principle of the method will be described using the example where a foreground object partially occludes the background as discussed in the Background of the Invention. Referring now to FIG. 4, a digital image processing system 40 suitable for practicing the invention is shown. The system includes an image processing computer 26, such as a Silicon Graphics Indiga II or Sun Sparc Workstation a video display 24, keyboard 30, and an input device, such as a mouse 28, that is connected to the processing computer 26 to provide a graphical user interface for processing motion image sequences (frames of motion images). A digital image storage device 32, such as a magnetic hard drive or optical disc memory, stores the original motion image sequences and stores the processed motion image sequence(s).

Referring now to FIG. 5A, to process a motion image sequence according to the present invention, the first step 42 is to identify the foreground and the background objects in the First and Second image frames at times t.sub.P and t.sub.Q by interactively drawing the boundaries of the foreground object(s) using any well known graphic user interface that provides for outlining areas in an image. Next, in step 44 and 46 the computer 26 estimates the velocity vector fields for the foreground and the background, respectively, in the Intermediate Frame at time t. The computer 26, then in step 48, identifies, segments, and labels the foreground and background regions for the intermediate image of the Intermediate Frame in the following manner: For each pixel determination in the intermediate image, the computer 26 uses the velocity vector at the pixel's location to locate corresponding points in Frames 1 and 2, at times t.sub.P and t.sub.Q. If the pixel's location in both frames is in the foreground, the pixel is identified as a foreground pixel and labeled a 2, per step 50, otherwise, the pixel is identified as a background pixel and is labeled as a 1, per step 52. With overlapping foreground pixels, the first overlapping pixels are labeled 3 with overlaps on overlaps the label number increases by 1 for each overlap. Referring now to FIG. 5B, in step 54 the computer 26 obtains the intermediate image for time t as follows: For each pixel in the intermediate image, if it is a foreground pixel, the computer uses the estimated foreground velocity vector at that pixel to find the correspondingly located pixel in the two frames, and computes the time weighted average of the two points. If the pixel is a background pixel, the computer finds the pixel's value in the two frames using the estimated background velocity vector and if they are both background pixel values, computes the time weighted average of the two points. If only one of the two corresponding pixel values is a background pixel value, the computer 26 uses that value as the background pixel value; and if they are both foreground the computer determines the value of the pixel from the pixel values of neighboring background pixels in the intermediate image. For purposes of illustration, the Intermediate image Frame in FIGS. 2 and 3 has been shown at a time t midway between Frame 1 and Frame 2. It should be noted that the method of the present invention may be employed to generate an intermediate image or a number of intermediate images at any time or times between any two captured frames.

Each pixel at the output of step 54 has one assigned estimated velocity vector field value. Those pixels are next processed in step 56 to interpolate the values of each pixels' neighbors using the label, the estimated velocity vector field value, and the pixels' image value.

The output of step 56 is a set of interpolated Intermediate Frame pixel values 58 that may be used in conjunction with an output device 60, such as a display device, or stored in a storage device, or transmitted to a remote location for viewing, modification, storing, or printing.

To simplify and in turn speed-up the computer's processing time, it is advantageous to only process the pixel's illuminance values and to set the R(red) and G(green) pixel values to some constant, for example, zero for all pixels outside of the encompassed objects of interest. In other words the image outside of objects having motion from one frame to another will appear blue, if the R and G pixel values are set to zero or to some other uniform color if they are all set to other values.

Referring now to FIG. 6A, the preferred embodiment of the method of the present invention will now be described in more detail. Using a graphical user interface, the user associates pixels of the existing frames with a "label" and a "flag." The (missing) pixels of the Intermediate Frame are associated with a "label" only. Labels are given integer values 1, 2, 3, . . . depending on the number of overlapping objects appearing in a frame. Flags represent either a "valid" or an "invalid" value, describing the validity of optical flow constraint as will be clear in the discussion that follows.

Referring to FIG. 6A, the truck, object 10, is encompassed (outlined) using any of the well known programs previously referenced to form the outline 70. The non-moving background to be occluded is the road portion labeled 16, as the truck 10 moves in the direction labeled D. Viewing FIG. 6A in conjunction with FIG. 6B, the occluding foreground (the truck) is assigned a label "2" and the background is assigned a label "1." The convention is such that the occluding object (in this case 10 and 10') is assigned a larger label than the occluded region. The user identifies the occluded regions on Frames 1 and 2 at times t.sub.P and t.sub.Q, respectively. For example, the objects 10 and 10' are assigned the label "2" and the background 12 is assigned a label "1." Determination of the corresponding labels in the Intermediate Frame is discussed below.

This implementation is described generally as having the following steps:

Step 1. Create two frames from the original frames, at time t.sub.P and t.sub.Q, by setting the image intensity to zero for regions labeled as "1." (Actual image intensity is maintained at pixels labeled as "2".) This effectively separates the foreground object from the background. The edges of the foreground object are now clearly defined for the Fogel application of the algorithm, without any interference from the background detail. The flags corresponding to the pixels of these two frames are then set to "valid" as previously states, with a color image, for example, an RGB image, the pixel values representing red and green may be set to zero and the blue pixel values can be set to a constant blue value. With that approach only the pixel values within the outlined object will be in R, G, B color and the background will be a substantially constant blue color.

Step 2. At this step, motion of the foreground object (truck) is estimated as described in detail below: estimate the linear motion trajectories between the frames created at the previous step, passing through a predetermined subset of pixel locations at t, using the Fogel algorithm of Appendix A. The resulting motion vector field is designated by .upsilon..sub.2.sup.t.

Step 3. Create two frames from the original frames, at times t.sub.P and t.sub.Q, by setting the flags for pixels corresponding to label "2" to "invalid." Flags for the reset of the pixels are set to "valid."

Step 4. At this step, motion of the background object is estimated as described in detail below: estimate the linear motion trajectories between the two frames created at the previous step, passing through the same pixel grid at t, using the Fogel algorithm. During estimation, suppress the optical flow constraint completely unless the trajectory intercepts regions with flags set to "invalid". The resulting motion vector field is denoted by .upsilon..sub.1.sup.t. At the end of this step, at every pixel grid location at t, there are two trajectories/vectors .upsilon..sub.1.sup.t and .upsilon..sub.2.sup.t.

Step 5. The estimated vector fields .upsilon..sub.1.sup.t and .upsilon..sub.2.sup.t both are interpolated to form dense vector fields where every pixel location, (i,j) at t has two associated vectors, .upsilon..sub.1.sup.t (i,j) and .upsilon..sub.2.sup.t (i,j).

Step 6. At this step, the intermediate frame at time t is labeled: for every pixel grid location at t, if the associated trajectory intercepts regions labeled as "2" in both existing frames, we label that pixel as "2". Otherwise it is labeled as "1." An ideal labeling of the intermediate frame is shown in FIG. 6B. According to this labeling, each pixel is associated with a motion vector (or trajectory); pixel at location (i,j) is associated with vector .upsilon..sub.1.sup.t (i,j) if it is labeled as "1," otherwise .upsilon..sub.2.sup.5 (i,j) is the trajectory passing through (i,j). For every grid location, the appropriate vector is maintained.

Step 7. In this step, the covered and uncovered regions at the intermediate frame are detected and handled by appropriate interpolation. We denote the labels of co-located pixels at positions (i,j), at times t.sub.P, t and t.sub.Q, as L0(i,j).rarw.L(i,j).fwdarw.L1(i,j). We have the following rules for appropriate interpolation:

IF 2.rarw.2.fwdarw.2 THEN use bilinear interpolation of the corresponding image values at times t.sub.P and t.sub.Q, on the trajectory .upsilon..sub.2.sup.t (i,j) passing through (i,j), to fill in the pixel value at location (i,j) at time t. The weights in the bilinear interpolation depend on the relative time-distance between the intermediate and the existing frames. The trajectory may point to interpixel values at the existing frames. Image vectors at such locations are computed via spatial bilinear interpolation of neighboring pixel values.

IF 1.rarw.1.fwdarw.1 THEN use bilinear interpolation of the corresponding image values at times t.sub.P and t.sub.Q, on the trajectory .upsilon..sub.1.sup.t (i,j) passing through (i,j) to fill in the pixel value at location (i,j) at time t.

IF 2.rarw.1.fwdarw.1 THEN pixel location (i,j) is within an uncovered region at the intermediate frame. (This can be seen in FIG. 6). This region should therefore be a part of the background. The image intensity at (i,j) should be taken from the future frame at t.sub.Q, since this pixel is covered in the previous frame at time t.sub.P. The intensity at (i,j) at time t is obtained by copying the intensity of the location where .upsilon..sub.1.sup.t (i,j) intercepts the frame at time t.sub.Q. Again, if the intercept is positioned at an inter-pixel location, the intensity value is obtained using spatial bilinear interpolation.

IF 1.rarw.1.fwdarw.2 THEN pixel location (i,j) is within a covered region at the intermediate frame. (This can be seen in FIG. 6). This region should therefore be a part of the background. However, the image intensity at (i,j) should be taken from the previous frame at t.sub.P, since this pixel is covered in the future frame at time t.sub.Q. The intensity at (i,j) at time t.sub.k is obtained by copying the intensity of the location where .upsilon..sub.1.sup.t (i,j) intercepts the frame at time t.sub.P.

In order to reduce computation in Step 2 above, vector field .upsilon..sub.2.sup.t can be estimated only for those pixel grid positions that lie in the convex hull of regions labeled as "2" in frames at time t.sub.P and t.sub.Q.

The above method can be extended to an arbitrary number of occluding objects. Suppose that there are three objects, where object "3" (labeled as "3") is partially occluding object "2," and object "2" is partially occluding the background object labeled as "1." In this case, three motion vector fields, .upsilon..sub.1.sup.t, .upsilon..sub.2.sup.t, .upsilon..sub.3.sup.t, are estimated at the intermediate frame. In estimating .upsilon..sub.1.sup.t, pixels with labels smaller than "3" are set to zero, and all flags are set to "valid". For .upsilon..sub.2.sup.t, pixels with labels smaller than "2" (i.e. "1") are set to zero, the flags for pixels with larger labels (i.e. "3") are set to "invalid."

The object based approach to detection and handling of occlusions in motion-compensated frame interpolation according to the present invention requires the user to input boundaries of occluding objects through a graphical interface. It is useful, for processing difficult cases, to create smooth slow and fast motion. The method is best used by applying it locally to difficult image areas only, rather than to the entire frames.

An even more detailed description of the embodiment outlined above, sufficient to enable a person of ordinary skill in the art to program a digital computer to practice the method of the present invention is contained in the following mathematical based description:

The invention has been described with reference to a preferred embodiment. However, it will be appreciated that variations and modifications can be effected by a person of ordinary skill in the art without departing from the scope of the invention.

We shall deal with images of a scene taken at time moments within some time interval t. 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. Given a time-varying image sequence of a scene corresponding to some sequence of time moments from the time interval t, it is often desirable to construct a time-varying image sequence of the same scene corresponding to some other sequence of time moments from the time interval t by constructing an intermediate image(s) (the images of a desired time-varying image sequence not present in a given time-varying image sequence). The algorithm to be described relates to a method for constructing a desired time-varying image sequence from a given time-varying image sequence by constructing an intermediate image herein after called the missing image.

To construct a missing image, first the correspondence between the missing image and the images of a given time-varying image sequence that surround the missing image is established. Then this correspondence is used to estimate the missing image through an interpolation process. The correspondence between the missing image and the surrounding images of a given time-varying image sequence establishes the apparent movements in these three images caused by the movements of objects in the scene. These apparent movements can be characterized by a velocity vector field defined on the missing image and linearly relating each point of the missing image to the points of the surrounding images. The construction of the missing image is then done in two stages: first the velocity vector field is estimated, second the estimate of the velocity vector field is used in a motion-compensated interpolation to obtain the missing image. Not all image points (pixels) will have corresponding points (pixels) in the given surrounding images, primarily due to occlusions, and allowances must be made to construct missing images with these occurrences.

Every image of the given time-varying image sequence is assumed to be obtained by projecting the light reflected by the visible portions of the objects in the scene at the time the image is taken onto an image plane, and then identifying the irradiance values of each image point (pixel) in the image plane. 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 given time-varying image sequence can be used to estimate the velocity vector fields of the images. The original scene is shown captured in Frame 1 as a plurality of pixel values, P.sub.1 -P.sub.3, for example. As is well known, the number of pixels to represent an image is dependent on a number of factor including the number of persons used in the photo sensor array capturing the sequence of scenes. In a three color system, three pixel values are used to record the scenes color at a pixel paint (location). These are denoted as the red, green, and blue values for an R, G, B, system.

An image of a scene is modeled as a perspective projection of a scene onto an image plane and captured as a frame. The axes of the scene coordinate system are defined as the horizontal axis, the vertical axis, and the depth axis. The images of the scene are assumed to be in the same plane which is parallel to the horizontal and the vertical axes of the scene coordinate system. It is assumed that the horizontal and the vertical axes of each image are parallel respectively to the horizontal and the vertical axes of the scene. In what follows the horizontal image coordinates are denoted with the symbol x, the vertical image coordinates are denoted with the symbol y, and the time moments are denoted with the symbol t. For the missing image at the time moment t the horizontal component of the estimate of the velocity vector at the image point (x,y) is denoted by u(x,y,t), while the vertical component of the estimate of the velocity vector at the image point (x,y) is denoted by .upsilon.(x,y,t). The estimate of the velocity vector field corresponding to the missing image at the time t is defined on some grid of points .OMEGA..sub.t in the image plane and is denoted as (u(t), v(t)).ident.{(u(x,y,t), .upsilon.(x,y,t)).vertline.(x,y).epsilon..OMEGA..sub.t }. The grid of image points .OMEGA..sub.t is assumed to be generated by some vectors h.sub.t,1 and h.sub.t,2 in the image plane, which means that every image point (x,y) from the grid .OMEGA..sub.t can be expressed in the form (x,y)=ih.sub.t,1 +jh.sub.t,2, where ih.sub.t,1 +jh.sub.t,2 is a linear combination of the vectors h.sub.t,1, h.sub.t,2 with integer coefficients i,j.

Let the time moments of the given time-varying image sequence be t.sub.P, t.sub.Q ; and let the time moments of the desired time-varying image sequence be t.sub.1, . . . , t.sub.K. The sequence of time moments t.sub.1, . . . t.sub.K is assumed to be increasing and intermediate to the increasing sequence of time moments t.sub.P, t.sub.Q. The given image of the scene taken at the time moment t.sub.P, called the previous image, and the given image of the scene taken at the time moment t.sub.Q, called the next image, are digitized into multi-banded digital images, RGB for example, with three bands, B. For a monochrome version there would only be one band B. Each band of the previous and next digital images is defined as a two-dimensioned array of numbers representing for example, the digital values of R, G and B, known as pixel values, taken from a grid of points on the image plane. The bands of the previous digital image are defined on some grid of image points .OMEGA..sub.t.sbsb.P, which is assumed to be generated by some vectors h.sub.t.sbsb.P,.sub.1 and h.sub.t.sbsb.P,.sub.2 in the image plane. The bands of the next digital image are defined on some grid of image points .OMEGA..sub.t.sbsb.Q, which is assumed to be generated by some vectors h.sub.t.sbsb.Q,.sub.1 and h.sub.t.sbsb.Q,.sub.2 in the image plane.

Given the previous digital image, which is defined on the grid of image points .OMEGA..sub.t.sbsb.P, taken at the time moment t.sub.P and consisting of the bands I.sub.t.sbsb.P,.sub.b, b=1, . . . , B, where b, for example, at value 1 represents the red pixel component value, and b at 2 the green component value, and b at 3 the blue component value, and given the next digital image, which is defined on the grid of image points .OMEGA..sub.t.sbsb.Q, taken at the time moment t.sub.Q and consisting of the bands I.sub.t.sbsb.Q,.sub.b, b=1, . . . , B, the process of constructing the desired time-varying image sequence consisting of the bands I.sub.t.sbsb.1,.sub.b, . . . , I.sub.t.sbsb.K,.sub.b, b=1, . . . , B can be described as follows. First, the estimates (u(t.sub.1),v(t.sub.1)), . . . , (u(t.sub.K),v(t.sub.K)) of the velocity vector fields, corresponding to the missing images at the time moments t.sub.1, . . . , t.sub.K respectively, are computed. The estimates (u(t.sub.1),v(t.sub.1)), . . . , (u(t.sub.K),v(t.sub.K)) of the velocity vector fields are defined on the grids of image points .OMEGA..sub.t.sbsb.1, . . . , .OMEGA..sub.t.sbsb.K respectively. Then for each band b=1, . . . , B and for each time moment t.sub.k, k=1, . . . , K, the estimate I.sub.t.sbsb.k,.sub.b of the missing image band, which is defined on the grid of image points .OMEGA..sub.t.sbsb.k, is obtained by the motion-compensated interpolation process.

It is an object of the invention to improve the accuracy of the estimation of the velocity vector field (u(t.sub.1),v(t.sub.1)), . . . , (u(t.sub.K),v(t.sub.K)), by the following actions: first, partitioning the grids .OMEGA..sub.t.sbsb.P, .OMEGA..sub.t.sbsb.Q into corresponding disjoint regions, .OMEGA..sub.t.sbsb.P.sub..function., .OMEGA..sub.t.sbsb.Q.sub..function., .function.=1, . . . , F, representing distinct objects in the images of a scene; second, performing the estimation process separately for each pair .OMEGA..sub.t.sbsb.P.sub..function., .OMEGA..sub.t.sbsb.Q.sub..function. of the corresponding regions, resulting in the velocity vector fields (u.sub..function. (t.sub.1),v.sub..function. (t.sub.1)), . . . , (u.sub..function. (t.sub.K),v.sub..function. (t.sub.K)), .function.=1, . . . , F; and then combining the velocity vector fields (u.sub.1 (t.sub.k),v.sub.1 (t.sub.k)), . . . , (u.sub.F (t.sub.k),v.sub.F (t.sub.k)) to form the velocity vector field (u(t.sub.k),v(t.sub.k)), for every k=1, . . . , K. The partitioning is specified by providing the corresponding sets of connected line segments on the previous and the next images representing the boundaries of the regions. This task is performed by an operator outlining the objects using a graphical user interface to encompass each object. The partitioning is accomplished by integrating the specified boundaries of the previous and next images into the regions and then attaching a distinct integer label to the points of each set of the corresponding regions. The following convention is used: the regions with the higher label values may occlude the regions with the lower label values, and the regions with the higher label values may not be occluded by the regions with the lower label values.

The proposed method of estimating the velocity vector field is based on the multi-level resolution approach where for every label the estimation process starts at the coarsest level of resolution with some constant velocity vector field taken as the initial estimate. This initial estimate of the velocity vector field is iteratively improved by solving the system of nonlinear equations and the resulted improved estimate is projected to the next finer level of resolution. The process is repeated until the final finest resolution level is reached.

In what follows, the symbol .sigma. is used to denote the current resolution level, the symbol M is used to denote the number of resolution levels, and the sequence of symbols .sigma..sub.1, . . . , .sigma..sub.M is used to denote the sequence of resolution levels listed in the fine-to-coarse order. For each label .function.=1, . . . , F, each resolution level .sigma.=.sigma..sub.1, . . . , .sigma..sub.M, and each k=1, . . . , K, the estimate of the velocity vector field at the resolution level .sigma. corresponding to the label .function. and the missing image at the time moment t.sub.k is defined on some grid .OMEGA..sub.t.sbsb.k.sup..sigma. of points in the image plane and is denoted as (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)).ident.{(u.sub.f.sup..sigma. (x,y,t.sub.k),.upsilon..sub..function..sup..sigma. (x,y,t.sub.k)).vertline.(x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma. }, while the improvement of the estimate of the velocity vector field at the resolution level .sigma. corresponding to the label .function. and the missing image at the time moment t.sub.k is defined on the same grid .OMEGA..sub.t.sbsb.k.sup..sigma. of points in the image plane and is denoted as (.DELTA.u.sub.f.sup..sigma. (t.sub.k),.DELTA.v.sub.f.sup..sigma. (t.sub.k)).ident.{(.DELTA.u.sub.f.sup..sigma. (x,y,t.sub.k),.DELTA..upsilon..sub.f.sup..sigma. (x,y,t.sub.k)).vertline.(x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma. }. For each resolution level .sigma.=.sigma..sub.1, . . . , .sigma..sub.M, and each k=1, . . . , K, the following is assumed to be held; the grid of points .OMEGA..sub.t.sbsb.k,.sup..sigma. is generated by some vectors h.sub.t.sbsb.k.sup..sigma..sub.1 and h.sub.t.sbsb.k,.sup..sigma..sub.2 in the image plane; and the grid of points .OMEGA..sub.t.sbsb.k in the image plane is assumed to be surrounded by the grid of points .OMEGA..sub.t.sbsb.k.sup.94 , meaning that the convex hall of the grid of points .OMEGA..sub.t.sbsb.k belongs to the convex hall of the grid of points .OMEGA..sub.t.sbsb.k.sup..sigma..

For each label f=1, . . . , F, on each level .sigma.=.sigma..sub.1, . . . , .sigma..sub.M of the multi-level resolution process, and for each k=1, . . . , K, a system of nonlinear equations relative to the unknown current estimate (u.sub..function..sup..sigma. (t.sub.k), v.sub..function..sup..sigma. (t.sub.k)) of the velocity vector field is formed, then the iterative improvement (.DELTA.u.sub.f.sup..sigma. (t.sub.k),.DELTA.v.sub.f.sup..sigma. (t.sub.k)) of the current estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field is achieved by solving the system of linear equations which is the linearization of this system of nonlinear equations about the current estimate (u.sub.f.sup..sigma. (t.sub.k),(v.sub.f.sup..sigma. (t.sub.k) of the velocity vector field. The iterative improvement (.DELTA.u.sub.f.sup..sigma. (t.sub.k),.DELTA.v.sub.f.sup..sigma. (t.sub.k) of the current estimate (u.sub.f.sup..sigma. (t.sub.k),(v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field is added to the current estimate (u.sub.f.sup..sigma. (t.sub.k),(v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field to obtain a new current estimate of the velocity vector field and the above process of the iterative improvement of the current estimate (u.sub.f.sup..sigma. (t.sub.k),(v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field is repeated until the system of nonlinear equations relative to the current estimate (u.sub.f.sup..sigma. (t.sub.k),(v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field is properly satisfied at which point the current estimate (u.sub.f.sup..sigma. (t.sub.k),(v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field is taken as the final estimate (u.sub.f,.sup..sigma..sub.1 (t.sub.k),(v.sub.f,.sup..sigma..sub.1 (t.sub.k)) of the velocity vector field at the resolution level .sigma..

For every f=1, . . . , F, every .sigma.=.sigma..sub.1, . . . , .sigma..sub.M, and every k=1, . . . , K, the system of nonlinear equations relative to the unknown current estimate (u.sub.f.sup..sigma. (t.sub.k),(v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field is defined in terms of the unknown current estimate (u.sub.f.sup..sigma. (t.sub.k),(v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field, and in terms of the quantities related to the values of previous and next images of the given time-varying image sequence and to the values of the generalized spatial partial derivatives of previous and next images of the given time-varying image sequence. The quantities related to the values of the previous image of the given time-varying image sequence, called the previous filtered image function and denoted as g.sup..sigma. (x,y,t.sub.P), and the quantities related to the values of the next image of the given time-varying image sequence, called the next filtered image function and denoted as g.sup..sigma. (x,y,t.sub.Q), are obtained by selectively, with respect to the partitioning, filtering each band of previous and next images with two filters: the horizontal derivative of the Gaussian-like filter and the vertical derivative of the Gaussian-like filter. The selective, with respect to the partitioning, filtering of previous and next images is defined in the next section as the filtering that does not mix the image values at the image points belonging to different regions of the respective grid partitioning. The quantities related to the values of the generalized x-spatial partial derivative of the previous image of the given time-varying image sequence, called the x-spatial partial derivative of the previous filtered image function and denoted as g.sub.x.sup..sigma. (x,y,t.sub.P), and the quantities related to the values of the generalized x-spatial partial derivative of the next image of the given time-varying image sequence, called the x-spatial partial derivative of the next filtered image function and denoted as g.sub.x.sup..sigma. (x,y,t.sub.Q), are obtained by selectively, with respect to the partitioning, filtering each band of previous and next images with the x-spatial partial derivative of each of the above two filters. The quantities related to the values of the generalized y-spatial partial derivative of the previous image of the given time-varying image sequence, called the y-spatial partial derivative of the previous filtered image function and denoted as g.sub.y.sup..sigma. (x,y,t.sub.P), and the quantities related to the values of the generalized y-spatial partial derivative of the next image of the given time-varying image sequence, called the y-spatial partial derivative of the next filtered image function and denoted as g.sub.y.sup..sigma. (x,y,t.sub.Q), are obtained by selectively, with respect to the partitioning, filtering each band of previous and next images with the y-spatial partial derivative of each of the above two filters.

For each t=t.sub.P, t.sub.Q, the grid of all points on the image plane generated by the vectors h.sub.t,1, h.sub.t,2 is denoted by the symbol H(h.sub.t,1,h.sub.t,2). The symbol G is used for the set of indices with each index g.epsilon.G specifying a particular band and a particular filter. For each time moment t=t.sub.P, t.sub.Q, each resolution level .sigma.=.sigma..sub.1, . . . , .sigma..sub.M, and each index g.epsilon.G, the filter function is denoted by the symbol .phi..sub.g.sup..sigma. (x,y,t). In the case of the horizontal derivative of the Gaussian-like filter, the filter function .phi..sub.g.sup..sigma. (x,y,t) is defined by the relation ##EQU1## for every (x,y).epsilon.H(h.sub.t,1 h.sub.t,2), such that .alpha..sub..sigma..sup.2 x.sup.2 +.beta..sub..sigma..sup.2 y.sup.2 <.omega..sup.2 ; while in the case of the vertical derivative of the Gaussian-like filter, the filter function .phi..sub.g.sup..sigma. (x,y,t) is defined by the relation ##EQU2## for every (x,y).epsilon.H(h.sub.t,1,h.sub.t,2), such that .alpha..sub..sigma..sup.2 x.sup.2 +.beta..sub..sigma..sup.2 y.sup.2 <.omega..sup.2. Here .alpha..sub..sigma., .beta..sub..sigma. and .omega. are positive constants. At the coarsest resolution level the constants .alpha..sub..sigma., .beta..sub..sigma. are set to a relatively small value, and at each finer resolution level the constants .alpha..sub..sigma., .beta..sub..sigma. are set to the increased by a fixed multiplicative factor values of the constants .alpha..sub..sigma., .beta..sub..sigma. at the coarser resolution level.

For each .sigma.=.sigma..sub.1, . . . , .sigma..sub.M, the previous filtered image functions g.sup..sigma. (x,y,t.sub.P), g.epsilon.G, the x-spatial partial derivatives of the previous filtered image functions g.sub.x.sup..sigma. (x,y,t.sub.P), g.epsilon.G, and the y-spatial partial derivatives of the previous filtered image functions g.sub.y.sup..sigma. (x,y,t.sub.P), g.epsilon.G are defined on the subgrid .OMEGA..sub.t.sbsb.P.sup..sigma. of the grid of points .OMEGA..sub.t.sbsb.P on the image plane, which is generated by some vectors h.sub.t.sbsb.P,.sup..sigma..sub.1 and h.sub.t.sbsb.P,.sup..sigma..sub.2, with the property: for every (x.sub.1,y.sub.1).epsilon..OMEGA..sub.t.sbsb.P.sup..sigma. and every (x.sub.2,y.sub.2).epsilon.H(h.sub.t.sbsb.P,.sup..sigma..sub.1,h.sub.t.sbsb .P,.sup..sigma..sub.2), such that .alpha..sub..sigma..sup.2 x.sub.2.sup.2 +.beta..sub..sigma..sup.2 y.sub.2.sup.2 <.omega..sup.2, the point (x,y)=(x.sub.1,y.sub.1)+(x.sub.2,y.sub.2) belongs to the grid .OMEGA..sub.t.sbsb.P. For each .sigma.=.sigma..sub.1, . . . , .sigma..sub.M, the next filtered image functions g.sup..sigma. (x,y,t.sub.Q), g.epsilon.G, the x-spatial partial derivatives of the next filtered image functions g.sub.x.sup..sigma. (x,y,t.sub.Q), g.epsilon.G, and the y-spatial partial derivatives of the next filtered image functions g.sub.y.sup..sigma. (x,y,t.sub.Q), g.epsilon.G are defined on the subgrid .OMEGA..sub.t.sbsb.Q.sup..sigma. of the grid of points .OMEGA..sub.t.sbsb.Q on the image plane, which is generated by some vectors h.sub.t.sbsb.Q,.sup..sigma..sub.1 and h.sub.t.sbsb.Q,.sup..sigma..sub.2, with the property: for every (x.sub.1,y.sub.1).epsilon..OMEGA..sub.t.sbsb.Q.sup..sigma. and every (x.sub.2,y.sub.2).epsilon.H(h.sub.t.sbsb.Q,.sup..sigma..sub.1,h.sub.t.sbsb .Q,.sup..sigma..sub.2), such that .alpha..sub..sigma..sup.2 x.sub.2.sup.2 +.beta..sub..sigma..sup.2 y.sub.2.sup.2 <.omega..sup.2, the point (x,y)=(x.sub.1,y.sub.1)+(x.sub.2,y.sub.2) belongs to the grid .OMEGA..sub.t.sbsb.Q.

The multi-level resolution process is build around a multi-level resolution pyramid meaning that on each coarser resolution level the previous filtered image functions and their spatial partial derivatives, the next filtered image functions and their spatial partial derivatives, and the estimates of the velocity vector fields are defined on subgrids of the corresponding grids at a finer resolution level. The grid rotation method, which complies with the requirements of the resolution pyramid, is used to define the grids .OMEGA..sub.t.sbsb.P.sup..sigma., .OMEGA..sub.t.sbsb.Q.sup..sigma., .OMEGA..sub.t.sbsb.k.sup..sigma., k=1, . . . , K, .sigma.=.sigma..sub.1, . . . , .sigma..sub.M. In this method, the grids are constructed sequentially starting with the finest resolution level and preceding towards the coarsest resolution level. The grids of the coarser resolution level are obtained from the grids of the finer resolution level by omiting the even grid points on the odd grid lines and by omiting the odd grid points on the even grid lines, and thereby rotating the grid lines of the grids of the finer resolution level by 45 degree. More precisely, these grids are defined as follows: for every t=t.sub.P, t.sub.Q, t.sub.I, . . . , t.sub.K, m=2, . . . , M, the grid .OMEGA..sub.t.sup..sigma..sbsp.m is sequentially defined as the subgrid of the grid .OMEGA..sub.t.sup..sigma..sbsp.m-1 generated by the vectors h.sub.t,1.sup..sigma..sbsp.m, h.sub.t,2.sup..sigma..sbsp.m, where the vectors h.sub.t,1.sup..sigma..sbsp.m, h.sub.t,2.sup..sigma..sbsp.m are defined by the relation

The proposed method of computing the estimates of the velocity vector fields, (u(t.sub.k),v(t.sub.k)), k=1, . . . , K, can be summerized as the follows:

1. Start with the label value F taken as the current label value f.

2. Start with the coarsest resolution level .sigma..sub.M taken as the current resolution level .sigma..

3. In the case of the current resolution level .sigma. being the coarsest resolution level .sigma..sub.M, set the initial estimate (u.sub.f,0.sup..sigma. (t.sub.k),v.sub.f,0.sup..sigma. (t.sub.k)) of the velocity vector field to some given (for example, constant) value, for every k=1, . . . , K; otherwise project the final estimate (u.sub.f,1.sup..delta.(.sigma. (t.sub.k),v.sub.f,1.sup..delta.(.sigma.) (t.sub.k)) of the velocity vector field at the previous coarser resolution level .delta.(.sigma.)=.sigma..sub.m+1 to the current resolution level .sigma.=.sigma..sub.m to obtain the initial estimate (u.sub.f,0.sup..sigma. (t.sub.k),v.sub.f,0.sup..sigma. (t.sub.k)) of the velocity vector field at the current resolution level .sigma., for every k=1, . . . , K. Take the initial estimate (u.sub.f,0.sup..sigma. (t.sub.k),v.sub.f,0.sup..sigma. (t.sub.k)) of the velocity vector field at the current resolution level as the current estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field at the current resolution level .sigma., for every k=1, . . . , K.

4. Iteratively improve the current estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field to obtain the final estimate (u.sub.f,1.sup..sigma. (t.sub.k),v.sub.f,1.sup..sigma. (t.sub.k)) of the velocity vector field at the current resolution level .sigma., for every k=1, . . . , K.

5. Take the next finer resolution level .sigma..sub.m-1 as the current resolution level .sigma..

6. Repeat step 3 to step 5 until the finest resolution level .sigma.=.sigma..sub.1 is reached.

7. Take the next smaller label value f-1 as the current label value f.

8. Repeat step 2 to step 7 until the smallest label value f=1 is reached.

9. Bilinearly interpolate the final estimate (u.sub.f,1.sup..sigma. (t.sub.k),v.sub.f,1.sup..sigma. (t.sub.k)) of the finest resolution level .sigma.=.sigma..sub.1 to the grid of points .OMEGA..sub.t.sbsb.k on the image plane to obtain the estimate (u.sub.f (t.sub.k),v.sub.f (t.sub.k)) of the velocity vector field corresponding to the label value f and to the missing image at the time moment t.sub.k, for every f=1, . . . , F, and every k=1, . . . K.

10. Combine the estimates (u.sub.f (t.sub.k),v.sub.f (t.sub.k)), f=1, . . . , F of the velocity vector field to obtain the estimate (u(t.sub.k), v(t.sub.k)) of the velocity vector field, for every k=1, . . . , K.

For a given time moment t=t.sub.P, t.sub.Q, and for every resolution level .sigma.=.sigma..sub.1, . . . , .sigma..sub.M, let the sequence of the subgrids .OMEGA..sub.t,f.sup..sigma., f=1, . . . , F of the grid .OMEGA..sub.t.sup..sigma. be the disjoint partitioning of the grid .OMEGA..sub.t.sup..sigma. defined in terms of the disjoint partitioning .OMEGA..sub.t,f, f=1, . . . , F of the grid .OMEGA..sub.t as follows. For every f=1, . . . , F, let .GAMMA..sub.t,f.sup..sigma. be the subgrid of the grid .OMEGA..sub.t.sup..sigma. consisting of the image points (x,y)=(x.sub.1,y.sub.1)+(x.sub.2,y.sub.2), for every (x.sub.1,y.sub.1).epsilon..OMEGA..sub.t,f .andgate..OMEGA..sub.t.sup..sigma. and every (x.sub.2,y.sub.2).epsilon.H(h.sub.t,1.sup..sigma.,h.sub.t,2.sup..sigma.) such that .alpha..sub..sigma..sup.2 x.sub.2.sup.2 +.beta..sub..sigma..sup.2 y.sub.2.sup.2 <.omega..sup.2. Then for every f=1, . . . , F, .OMEGA..sub.t,f.sup..sigma. is defined as the subgrid of the grid .OMEGA..sub.t.sup..sigma. consisting of those points of the union of the subgrids .GAMMA..sub.t,f.sup..sigma., . . . , .GAMMA..sub.t,F.sup..sigma. that do not belong to the union of the subgrids .GAMMA..sub.t,f+1.sup..sigma., . . . , .GAMMA..sub.t,F.sup..sigma..

For each time moment t=t.sub.P, t.sub.Q, each resolution level .sigma.=.sigma..sub.1, . . . , .sigma..sub.M, each label f=1, . . . , F, each index g.epsilon.G, and each image point (x,y).epsilon..OMEGA..sub.t,f.sup..sigma., the value g.sup..sigma. (x,y,t) of the filtered image function is given by the relation ##EQU3## the value g.sub.x.sup..sigma. (x,y,t) of the x-spatial partial derivative of the filtered image function g.sup..sigma. (x,y,t) is given by the relation ##EQU4## the value g.sub.y.sup..sigma. (x,y,t) of the y-spatial partial derivative of the filtered image function g.sup..sigma. (x,y,t) is given by the relation ##EQU5## where: b.sub.g is the band corresponding to the index g, and the function I.sub.t,f,b.sbsb.g (x',y') is equal to the image function I.sub.t,b.sbsb.g (x',y') for every (x',y').epsilon..OMEGA..sub.t,f and is equal to zero otherwise.

For each k=1, . . . , K, each .sigma.=.sigma..sub.1, . . . , .sigma..sub.M, and each f=1, . . . , F, the system of nonlinear equations relative to the unknown current estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field forms a compromise between optical flow and directional smoothness constraints. The optical flow constraints, represented by optical flow functions and their spatial partial derivatives, relate the values of the next and previous filtered image functions at the corresponding image points. The directional smoothness constraints, represented by the directional smoothness functions, relate the values of neighboring velocity vectors.

For each index g.epsilon.G specifying a particular filter and a particular image band, for the following three functions defined on the grid of image points (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma. : an optical flow function g.sub.t.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)), a spatial partial with respect to the component u.sub.f.sup..sigma. (x,y,t.sub.k) of the velocity vector derivative g.sub.tu.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the optical flow function g.sub.t.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)), and a spatial partial with respect to the component v.sub.f.sup..sigma. (x,y,t.sub.k) of the velocity vector derivative g.sub.tv.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the optical flow function g.sub.t.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)). The value of the optical flow function g.sub.t.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) is equal to the scaled by the factor 1/(t.sub.Q -t.sub.P) difference between the estimated next filtered image function and the estimated previous filtered image function at the points (x,y) of the grid .OMEGA..sub.t.sbsb.k.sup..sigma. where both the estimate next filtered image function and the estimated previous filtered image function are defined and equal to zero otherwise. The value of the spatial partial with respect to the component u.sub.f.sup..sigma. (x,y,t.sub.k) of the velocity vector derivative g.sub.tu.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the optical flow function is equal to the sum of the scaled by the factor (t.sub.Q -t.sub.k)/(t.sub.Q -t.sub.P) spatial partial with respect to the component x derivative of the estimated next filtered image function and the scaled by the factor (t.sub.k -t.sub.P)/(t.sub.Q -t.sub.P) spatial partial with respect to the component x derivative of the estimated previous filtered image function at the points (x,y) of the grid .OMEGA..sub.t.sbsb.k.sup..sigma. where both the spatial partial with respect to the component x derivative of the estimated next filtered image function and the spatial partial with respect to the component x derivative of the estimated previous filtered image function are defined and equal to zero otherwise. The value of the spatial partial with respect to the component v.sub.f.sup..sigma. (x,y,t.sub.k) of the velocity vector derivative g.sub.tv.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the optical flow function is equal to the sum of the scaled by the factor (t.sub.Q -t.sub.k)/(t.sub.Q -t.sub.P) spatial partial with respect to the component y derivative of the estimated next filtered image function and the scaled by the factor (t.sub.k -t.sub.P)/(t.sub.Q -t.sub.P) spatial partial with respect to the component y derivative of the estimated previous filtered image function at the points (x,y) of the grid .OMEGA..sub.t.sbsb.k.sup..sigma. where both the spatial partial with respect to the component y derivative of the estimated next filtered image function and the spatial partial with respect to the component y derivative of the estimated previous filtered image function are defined and equal to zero otherwise.

To obtain the estimated next filtered image function, the spatial partial with respect to the component x derivative of the estimated next filtered image function, and the spatial partial with respect to the component y derivative of the estimated next filtered image function, find the grid of image points .OMEGA..sub.t.sbsb.Q.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)), which are the taken at the time moment t.sub.Q estimated perspective projections onto the image plane of the visible points in the scene, whose perspective projections onto the image plane at the time moment t.sub.k are the grid points .OMEGA..sub.t.sbsb.k.sup..sigma. ; and then bilinearly interpolate the next filtered image function g.sup..sigma. (x,y,t.sub.Q), the spatial partial with respect to the component x derivative g.sub.x.sup..sigma. (x,y,t.sub.Q) of the next filtered image function, and the spatial partial with respect to the component y derivative g.sub.y.sup..sigma. (x,y,t.sub.Q) of the next filtered image function, from the grid of points .OMEGA..sub.t.sbsb.Q.sup..sigma. to those points of the grid .OMEGA..sub.t.sbsb.Q.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) that are surrounded by the points of the grid .OMEGA..sub.t.sbsb.Q.sub.,f.sup..sigma.. To obtain the estimated previous filtered image function, the spatial partial with respect to the component x derivative of the estimated previous filtered image function, and the spatial partial with respect to the component y derivative of the estimated previous filtered image function, find the grid of image points .OMEGA..sub.t.sbsb.P.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)), which are the taken at the time moment t.sub.P estimated perspective projections onto the image plane of the visible points in the scene, whose perspective projections onto the image plane at the time moment t.sub.k are the grid points .OMEGA..sub.t.sbsb.k.sup..sigma. ; and then bilinearly interpolate the previous filtered image function g.sup..sigma. (x,y,t.sub.P), the spatial partial with respect to the component x derivative g.sub.x.sup..sigma. (x,y,t.sub.P) of the previous filtered image function, and the spatial partial with respect to the component y derivative g.sub.y.sup..sigma. (x,y,t.sub.P) of the previous filtered image function, from the grid of points .OMEGA..sub.t.sbsb.P.sup..sigma. to those points of the grid .OMEGA..sub.t.sbsb.P.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) that are surrounded by the points of the grid .OMEGA..sub.t.sbsb.P.sub.,f.sup..sigma.. To obtain the grid of image points .OMEGA..sub.t.sbsb.Q.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)), each image point of the grid .OMEGA..sub.t.sbsb.k.sup..sigma. is shifted by the amount equal to the scaled by the factor (t.sub.Q -t.sub.k) value of the estimated velocity vector defined at that grid point. To obtain the grid of image points .OMEGA..sub.t.sbsb.P.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)), each image point of the grid .OMEGA..sub.t.sbsb.k.sup..sigma. is shifted by the amount equal to the scaled by the factor (t.sub.P -t.sub.k) value of the estimated velocity vector defined at that grid point.

Let (x,y) be an image point from the grid .OMEGA..sub.t.sbsb.k.sup..sigma. where the estimated next filtered image function, its spatial partial with respect to the component x and to the component y derivatives, the estimated previous filtered image function, and its spatial partial with respect to the component x and to the component y derivatives are all defined, then the value g.sub.t.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k), v.sub.f.sup..sigma. (x,y,t.sub.k)) of the optical flow function, the value g.sub.tu.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the spatial partial with respect to the component u.sub.f.sup..sigma. (x,y,t.sub.k) of the velocity vector derivative of the optical flow function and the value g.sub.tv.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the spatial partial with respect to the component v.sub.f.sup..sigma. (x,y,t.sub.k) of the velocity derivative of the optical flow function at the image point (x,y) can be more precisely defined as follows. The point (x+(t.sub.Q -t.sub.k)u.sub.f.sup..sigma. (x,y,t.sub.k),y+(t.sub.Q -t.sub.k)v.sub.f.sup..sigma. (x,y,t.sub.k)) of the grid .OMEGA..sub.t.sbsb.Q.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) can be expressed in the following form

for some integer numbers i.sub.t.sbsb.Q.sub.,1.sup..sigma., i.sub.t.sbsb.Q.sub.,2.sup..sigma. and for some real numbers .theta..sub.t.sbsb.Q.sub.,1.sup..sigma., .theta..sub.t.sbsb.Q.sub.,2.sup..sigma. from the interval [0,1]. The fact that the estimated next filtered image function and its spatial partial with respect to the component x and to the component y derivatives are defined at the point (x,y) means that the point (x+(t.sub.Q -t.sub.k)u.sub.f.sup..sigma. (x,y,t.sub.k),y+(t.sub.Q -t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the grid .OMEGA..sub.t.sbsb.Q.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) is surrounded by the following points of the grid .OMEGA..sub.t.sbsb.Q.sub.,f.sup..sigma. :

The point (x+(t.sub.P -t.sub.k)u.sub.f.sup..sigma. (x,y,t.sub.k),y+(t.sub.P -t.sub.k)v.sub.f.sup..sigma. (x,y,t.sub.k)) of the grid .OMEGA..sub.t.sbsb.P.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) can be expressed in the following form

for some integer numbers i.sub.t.sbsb.P.sub.,1.sup..sigma., i.sub.t.sbsb.P.sub.,2.sup..sigma. and for some real numbers .theta..sub.t.sbsb.P.sub.,1.sup..sigma., .theta..sub.t.sbsb.P.sub.,2.sup..sigma. from the interval [0,1]. The fact that the estimated previous filtered image function and its spatial partial with respect to the component x and to the component y derivatives are defined at the point (x,y) means that the point (x+(t.sub.P -t.sub.k)u.sub.f.sup..sigma. (x,y,t.sub.k),y+(t.sub.P -t.sub.k)v.sub.f.sup..sigma. (x,y,t.sub.k)) of the grid .OMEGA..sub.t.sbsb.P.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) is surrounded by the following points of the grid .OMEGA..sub.t.sbsb.P.sub.,f.sup..sigma. :

The value g.sub.t.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the optical flow function at the image point (x,y) is given by the relation ##EQU6## where the value g.sup..sigma. (x+(t.sub.Q -t.sub.k)u.sub.f.sup..sigma. (x,y,t.sub.k),y+(t.sub.Q -t.sub.k).nu..sub.f.sup..sigma. (x,y,t.sub.k),t.sub.Q) of the estimated next filtered image function is given as the following bilinear interpolation ##EQU7## of the values of the next filtered image function at the image points defined by the relations (5-2), (5-3), (5-4), (5-5); while the value g.sup..sigma. (x+(t.sub.P -t.sub.k)u.sub.f.sup..sigma. (x,y,t.sub.k),y+(t.sub.P -t.sub.k)v.sub.f.sup..sigma. (x,y,t.sub.k),t.sub.P) of the estimated previous filtered image function is given as the following bilinear interpolation ##EQU8## of the values of the previous filtered image function at the image points defined by the relations (5-7), (5-8), (5-9), (5-10). The value g.sub.tu.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the spatial partial with respect to the component u.sub.f.sup..sigma. (x,y,t.sub.k) of the velocity vector derivative of the optical flow function at the image point (x,y) is given by the relation ##EQU9## where the value g.sub.x.sup..sigma. (x+(t.sub.Q -t.sub.k)u.sub.f.sup..sigma. (x,y,t.sub.k),y+(t.sub.Q -t.sub.k)v.sub.f.sup..sigma. (x,y,t.sub.k),t.sub.Q) of the spatial partial with respect to the component x derivative of the estimated next filtered image function is given as the following bilinear interpolation ##EQU10## of the values of the spatial partial with respect to the component x derivative of the next filtered image function at the image points defined by the relations (5-2), (5-3), (5-4), (5-5); while the value g.sub.x.sup..sigma. (x+(t.sub.P -t.sub.k)u.sub.f.sup..sigma. (x,y,t.sub.k),y+(t.sub.P -t.sub.k)v.sub.f.sup..sigma. (x,y,t.sub.k),t.sub.P) of the spatial partial with respect to the component x derivative of the estimated previous filtered image function is given as the following bilinear interpolation ##EQU11## of the values of the spatial partial with respect to the component x derivative of the previous filtered image function at the image points defined by the relations (5-7), (5-8), (5-9), (5-10). The value g.sub.tv.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the spatial partial with respect to the component v.sub.f.sup..sigma. (x,y,t.sub.k) of the velocity vector derivative of the optical flow function at the image point (x,y) is given by the relation ##EQU12## where the value g.sub.y.sup..sigma. (x+(t.sub.Q -t.sub.k)u.sub.f.sup..sigma. (x,y,t.sub.k),y+(t.sub.Q -t.sub.k)v.sub.f.sup..sigma. (x,y,t.sub.k),t.sub.Q) of the spatial partial with respect to the component y derivative of the estimated next filtered image function is given as the following bilinear interpolation ##EQU13## of the values of the spatial partial with respect to the component y derivative of the next filtered image function at the image points defined by the relations (5-2), (5-3), (5-4), (5-5); while the value g.sub.y.sup..sigma. (x+(t.sub.P -t.sub.k)u.sub.f.sup..sigma. (x,y,t.sub.k),y+(t.sub.P -t.sub.k)v.sub.f.sup..sigma. (x,y,t.sub.k),t.sub.P) of the spatial partial with respect to the component y derivative of the estimated previous filtered image function is given as the following bilinear interpolation ##EQU14## of the values of the spatial partial with respect to the components y derivative of the previous filtered image function at the image points defined by the relations (5-7), (5-8), (5-9), (5-10).

Each non-boundary grid point (x,y) of the grid .OMEGA..sub.t.sbsb.k.sup..sigma., k=1, . . . , K is surrounded by its eight nearest neighbors specifying eight different directions on the image plane. For each boundary grid point (x,y) of the grid .OMEGA..sub.t.sbsb.k.sup..sigma., k= 1, . . . , K the nearest neighbors are present only for some of these eight different directions. The symbol s is used to denote a vector on the image plane specifying one of these eight different directions on the image plane. The symbol S is used to denote the set of these eight vectors. Let .OMEGA..sub.t.sbsb.k.sub.,s.sup..sigma., k=1, . . . , K, s.epsilon.S be a subgrid of the grid .OMEGA..sub.t.sbsb.k.sup..sigma. with the property that every grid point (x,y) of the subgrid .OMEGA..sub.t.sbsb.k.sub.,s.sup..sigma. has the nearest in the direction s neighbor on the grid .OMEGA..sub.t.sbsb.k.sup..sigma.. For each k=1, . . . , K, and for each vector s from the set S specifying a particular direction on the image plane, for the directional smoothness function (s,.gradient.u.sub.f.sup..sigma. (x,y,t.sub.k)), (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma. for the horizontal component u.sub.f.sup..sigma. (x,y,t.sub.k) of the current estimate of the velocity vector field as the finite difference approximation to the directional derivative of the horizontal components u.sub.f.sup..sigma. (x,y,t.sub.k) of the current estimate of the velocity vector field as follows. For every point (x,y) of the subgrid .OMEGA..sub.t.sbsb.k.sub.,s.sup..sigma. of the grid .OMEGA..sub.t.sbsb.k.sup..sigma., the value of the directional smoothness function (s,.gradient.u.sub.f.sup..sigma. (x,y,t.sub.k)) at the point (x,y) is defined to be equal to the difference between the value u.sub.f.sup..sigma. (x.sub.k,s,y.sub.k,s,t.sub.k) of the horizontal component of the current estimate of the velocity vector at the nearest in the direction s neighbor (x.sub.k,s,y.sub.k,s) of this grid point (x,y) and the value u.sub.f.sup..sigma. (x,y,t.sub.k)) of the horizontal component of the current estimate of the velocity vector at this grid point (x,y). For every point (x,y) of the grid .OMEGA..sub.t.sbsb.k.sup..sigma. not belonging to the subgrid .OMEGA..sub.t.sbsb.k.sub.,s.sup..sigma. the value of the directional smoothness function (s,.gradient.u.sub.f.sup..sigma. (x,y,t.sub.k)) at the point (x,y) is defined to be equal to zero. For each k=1, . . . , K, and for each vector s from the set S specifying a particular direction on the image plane, form the directional smoothness function (s,.gradient.v.sub.f.sup..sigma. (x,y,t.sub.k)), (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma. for the vertical component .nu..sub.f.sup..sigma. (x,y,t.sub.k) of the current estimate of the velocity vector field as the finite difference approximations to the directional derivatives of the vertical component v.sub.f.sup..sigma. (x,y,t.sub.k) of the current estimate of the velocity vector field as follows. For every point (x,y) of the subgrid .OMEGA..sub.t.sbsb.k.sub.,s.sup..sigma. of the grid .OMEGA..sub.t.sbsb.k.sup..sigma., the value of the directional smoothness function (s,.gradient.v.sub.f.sup..sigma. (x,y,t.sub.k)) at the point (x,y) is defined to be equal to the difference between the value v.sub.f.sup..sigma. (x.sub.k,s,y.sub.k,s,t.sub.k) of the vertical component of the current estimate of the velocity vector at the nearest in the direction s neighbor (x.sub.k,s,,y.sub.k,s) of this grid point (x,y) and the value .nu..sub.f.sup..sigma. (x,y,t.sub.k)) of the vertical component of the current estimate of the velocity vector at this grid point (x,y). For every point (x,y) of the grid .OMEGA..sub.t.sbsb.k.sup..sigma. not belonging to the subgrid .OMEGA..sub.t.sbsb.k.sub.,s.sup..sigma. the value of the directional smoothness function (s,.gradient.v.sub.f.sup..sigma. (x,y,t.sub.k)) at the point (x,y) is defined to be equal to zero.

For every k=1, . . . , K, the system of nonlinear equations relative to the unknown current estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field is formed by combining the optical flow function g.sub.t.sup..sigma. and its spatial derivatives g.sub.tu.sup..sigma., g.sub.tv.sup..sigma. for each filter and each image band specified by the index g.epsilon.G together with the directional smoothness functions (s,.gradient.u.sub.f.sup..sigma.), (s,.gradient.v.sub.f.sup..sigma.) for each image direction s.epsilon.S and together with some constant parameters as follows ##EQU15## for every (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma.. Here: the expression q.sup.2 (.gradient.u.sub.f.sup..sigma.,.gradient.v.sub.f.sup..sigma.), (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma. represents a square of the norm of the gradient (.gradient.u.sub.f.sup..sigma.,.gradient.v.sub.f.sup..sigma.) of the velocity vector (u.sub.f.sup..sigma.,v.sub.f.sup..sigma.) and is defined by the relation ##EQU16## where .rho..sub.s is the length of the vector s, s.epsilon.S and p.sub.0, p.sub.1 are positive constant parameters; the expression b.sup.2 (s,.gradient.'g.sub.t.sup..sigma.), (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma., s.epsilon.S is defined by the relation ##EQU17## where

the expression (u.sub.f,0.sup..sigma.,v.sub.f,0.sup..sigma.).tbd.(u.sub.f,0.sup..sigma. (x,y,t.sub.k),v.sub.f,0.sup..sigma. (x,y,t.sub.k)), (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma. represents the initial estimate of the velocity vector at the image point (x,y); and r.sub.g, a.sub.s, c, b.sub.g, .gamma..sup..sigma. are positive constant parameters. The system of linear equations relative to the improvement (.DELTA.u.sub.f.sup..sigma. (t.sub.k),.DELTA.v.sub.f.sup..sigma. (t.sub.k)) of the current estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field is defined as the symmetric, positive definite, and diagonally dominant linearization of the system of nonlinear equations (6-1) about the current estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field as follows ##EQU18## for every (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma.. The solution of the system of linear equations (6-5) is obtained using a fast preconditioned iterative method which is based on the polynominal acceleration of the basic iterative method and can be described as follows.

For every k=1, . . . , K, and every f=1, . . . , F, the system of linear equations (6-5) can be expressed in the matrix notation as

where (.DELTA.u.sub.f.sup..sigma. (t.sub.k),.DELTA.v.sub.f.sup..sigma. (t.sub.k)).sup.T is the transpose vector of the row-vector (.DELTA.u.sub.f.sup..sigma. (t.sub.k),.DELTA.v.sub.f.sup..sigma. (t.sub.k)). The square matrix M.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) can be partitioned into four square submatrices corresponding to the vectors (.DELTA.u.sub.f.sup..sigma. (t.sub.k)) and (.DELTA.v.sub.f.sup..sigma. (t.sub.k)), respectively, so that the following relation is satisfied ##EQU19## Let the matrix D.sub.i,j.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) be the diagonal part of the matrix M.sub.i,j.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)), for every i,j=1,2, and let the matrix D.sup..sigma. (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)), be defined as in the relations ##EQU20## Then the preconditioned conjugate gradient iterative method for approximating the solution of the system of linear equations (7-1) can be described as follows.

For the sake of simplicity of the presentation the following notation shall be used:

Start with the initial approximation w.sub.0 which is identically equal to zero, and with the initial residual r.sub.0 which is identically equal to h, then for every n=0, 1, . . . , N-1 the following iterative procedure is applied:

The coefficient .alpha..sub.n is defined as in the relation ##EQU21## The vectors q.sub.n, p.sub.n are given by the relation

The coefficient .beta..sub.n and the vector z.sub.n are defined as in the relation ##EQU22##

For every time moment t.sub.k, k=1, . . . , K, let the sequence of the subgrids .OMEGA..sub.t.sbsb.k.sub.,f, f=1, . . . , F of the grid .OMEGA..sub.t.sbsb.k be the disjoint partitioning of the grid .OMEGA..sub.t.sbsb.k defined below in terms of the following: the disjoint partitioning .OMEGA..sub.t.sbsb.p.sub.,f, f=1, . . . , F of the grid .OMEGA..sub.t.sbsb.P, the disjoint partitioning .OMEGA..sub.t.sbsb.Q.sub.,f, f=1, . . . , F of the gird .OMEGA..sub.t.sbsb.Q, and the labeled estimate (u.sub.f (t.sub.k),v.sub.f (t.sub.k)), f=1, . . . , F of the velocity vector field. For every f=2, . . . , F, the subgrid .OMEGA..sub.t.sbsb.k.sub.,f of the grid .OMEGA..sub.t.sbsb.k consists of the image points (x,y) with the properties: i) for every point (x,y).epsilon..OMEGA..sub.t.sbsb.k.sub.,f, the point (x+(t.sub.P -t.sub.k)u.sub.f (x,y,t.sub.k),y+(t.sub.P -t.sub.k)v.sub.f (x,y,t.sub.k)) is not surrounded by the points of the subgrids .OMEGA..sub.t.sbsb.P.sub.,1, . . . , .OMEGA..sub.t.sbsb.P.sub.,f-1 ; ii) for every point (x,y).epsilon..OMEGA..sub.t.sbsb.k.sub.,f, the point (x+(t.sub.Q -t.sub.k)u.sub.f (x,y,t.sub.k),y+(t.sub.Q -t.sub.k)v.sub.f (x,y,t.sub.k)) is not surrounded by the points of the subgrids .OMEGA..sub.t.sbsb.Q.sub.,1, . . . , .OMEGA..sub.t.sbsb.Q.sub.,f-1 ; iii) in the case of the value of the label f being less than the value of the label F, every point (x,y) of the grid .OMEGA..sub.t.sbsb.k belonging to the subgrid .OMEGA..sub.t.sbsb.k.sub.,f does not belong to any of the subgrids .OMEGA..sub.t.sbsb.k.sub.,f+1, . . . , .OMEGA..sub.t.sbsb.k.sub.,F. The subgrid .OMEGA..sub.t.sbsb.k.sub.,1 of the grid .OMEGA..sub.t.sbsb.k consists of the image points (x,y) that do not belong to any of the subgrids .OMEGA..sub.t.sbsb.k.sub.,2, . . . , .OMEGA..sub.t.sbsb.k.sub.,F. Then for every label f=1, . . . , F and for every point (x,y).epsilon..OMEGA..sub.t,f the value of the estimate (u(x,y,t.sub.k),v(x,y,t.sub.k)) of the velocity vector field is defined to be equal to the value of the labeled estimate (u.sub.f (x,y,t.sub.k),v.sub.f (x,y,t.sub.k)) of the velocity vector field.

For each time moment t.sub.k, k=1, . . . , K, given the next digital image, which is defined on the grid of image points .OMEGA..sub.t.sbsb.Q, taken at the time moment t.sub.Q and consisting of the bands I.sub.t.sbsb.Q.sub.,1, . . . , I.sub.t.sbsb.Q.sub.,B, and given the previous digital image, which is defined on the grid of image points .OMEGA..sub.t.sbsb.P, taken at the time moment t.sub.P and consisting of the bands I.sub.t.sbsb.P.sub.,1, . . . , I.sub.t.sbsb.P.sub.,B, the method of estimating the missing digital image, which is defined on the grid of image points .OMEGA..sub.t.sbsb.k, taken at the time moment t.sub.k and consisting of the bands I.sub.t.sbsb.k.sub.,1, . . . , I.sub.t.sbsb.k.sub.,B based on the estimate (u(t.sub.k),v(t.sub.k)) of the velocity vector field, which is defined on the grid of image points .OMEGA..sub.t.sbsb.k, can be described as follows. For each image band b=1, . . . , B, the value of the estimate I.sub.t.sbsb.k.sub.,b of the missing digital image band is equal to the sum of the scaled by the factor (t.sub.k -t.sub.P)/(t.sub.Q -t.sub.P) value of the estimated next digital image band and the scaled by the factor (t.sub.Q -t.sub.k)/(t.sub.Q -t.sub.P) value of the estimated previous digital image band at the points (x,y) of the grid .OMEGA..sub.t.sbsb.k where both the estimated next digital image band and the estimated previous digital image band are defined. For each image band b=1, . . . , B, the value of the estimate I.sub.t.sbsb.k.sub.,b of the missing digital image band is equal to the value of the estimated next digital image band at the points (x,y) of the grid .OMEGA..sub.t.sbsb.k where only the estimated next digital image band is defined. For each image band b=1, . . . , B, the value of the estimate I.sub.t.sbsb.k.sub.,b of the missing digital image band is equal to the value of the estimated previous digital image band at the points (x,y) of the grid .OMEGA..sub.t.sbsb.k where only the estimated previous digital image band is defined. To obtain the estimated next digital image band, for every f=1, . . . , F find the grid of image points .OMEGA..sub.t.sbsb.Q.sub.,f (u(t.sub.k),v(t.sub.k)) which are the taken at the time moment t.sub.Q estimated perspective projections onto the image plane of the visible points in the scene, whose perspective projections onto the image plane at the time moment t.sub.k are the grid points .OMEGA..sub.t.sbsb.k.sub.,f ; and then bilinearly interpolate the belonging to the subgrid .OMEGA..sub.t.sbsb.k.sub.,f portion of the next digital image band I.sub.t.sbsb.Q.sub.,b from the grid of points .OMEGA..sub.t.sbsb.Q where the next digital image is defined to those points of the grid .OMEGA..sub.t.sbsb.Q.sub.,f (u(t.sub.k),v(t.sub.k)) that are surrounded by the points of the grid .OMEGA..sub.t.sbsb.Q.sub.,f. To obtain the estimated previous digital image band, for every f=1, . . . , F find the grid of image points .OMEGA..sub.t.sbsb.P.sub.,f (u(t.sub.k),v(t.sub.k)) which are the taken at the time moment t.sub.P estimated perspective projections onto the image plane of the visible points in the scene, whose perspective projections onto the image plane at the time moment t.sub.k are the grid points .OMEGA..sub.t.sbsb.k.sub.,f ; and then bilinearly interpolate the belonging to the subgrid .OMEGA..sub.t.sbsb.k.sub.,f portion of the previous digital image band I.sub.t.sbsb.P.sub.,b from the grid of points .OMEGA..sub.t.sbsb.P where the previous digital image is defined to those points of the grid .OMEGA..sub.t.sbsb.P.sub.,f (u(t.sub.k),v(t.sub.k)) that are surrounded by the points of the grid .OMEGA..sub.t.sbsb.P.sub.,f. To obtain the grid of image points .OMEGA..sub.t.sbsb.Q.sub.,f (u(t.sub.k),v(t.sub.k)), f=1, . . . , F, each image point of the grid .OMEGA..sub.t.sbsb.k.sub.,f is shifted by the amount equal to the scaled by the factor (t.sub.Q -t.sub.k) value of the estimated velocity vector defined at that grid point. To obtain the grid of image points .OMEGA..sub.t.sbsb.P.sub.,f (u(t.sub.k),v(t.sub.k)), f=1, . . . , F, each image point of the grid .OMEGA..sub.t.sbsb.k.sub.,f is shifted by the amount equal to the scaled by the factor (t.sub.k -t.sub.P) value of the estimated velocity vector defined at that grid point.

For a given label value f=1, . . . , F, let (x,y) be an image point from the grid .OMEGA..sub.t.sbsb.k.sub.,f, then for a given band value b=1, . . . , B the value of the missing digital image band I.sub.t.sbsb.k.sub.,b (x,y,u(x,y,t.sub.k),v(x,y,t.sub.k)) at the image point (x,y) can be more precisely defined as follows.

The point (x+(t.sub.Q -t.sub.k)u(x,y,t.sub.k),y+(t.sub.Q -t.sub.k)v(x,y,t.sub.k)) of the grid .OMEGA..sub.t.sbsb.Q.sub.,f (u(t.sub.k),v(t.sub.k)) can be expressed in the following form

for some integer numbers i.sub.t.sbsb.Q.sub.,1, i.sub.t.sbsb.Q.sub.,2 and for some real numbers .theta..sub.t.sbsb.Q.sub.,1, .theta..sub.t.sbsb.Q.sub.,2 from the interval [0,1]. If the estimated next digital image is defined at the point (x,y,), then the point (x+(t.sub.Q -t.sub.k)u(x,y,t.sub.k),y+(t.sub.Q -t.sub.k)v(x,y,t.sub.k)) of the grid .OMEGA..sub.t.sbsb.Q.sub.,f (u(t.sub.k),v(t.sub.k)) is surrounded by the following points of the grid .OMEGA..sub.t.sbsb.Q.sub.,f :

and the value I.sub.t.sbsb.Q.sub.,b (x+(t.sub.Q -t.sub.k)u(x,y,t.sub.k),y+(t.sub.Q -t.sub.k)v(x,y,t.sub.k)) of the estimated nest digital image band is given as the following bilinear interpolation ##EQU23## of the values of the next digital image band at the image points defined by the relations (9-2), (9-3), (9-4), (9-5).

The point (x+(t.sub.P -t.sub.k)u(x,y,t.sub.k),y+(t.sub.P -t.sub.k)v(x,y,t.sub.k)) of the grid .OMEGA..sub.t.sbsb.P.sub.,f (u(t.sub.k),v(t.sub.k)) can be expressed in the following form

for some integer numbers i.sub.t.sbsb.P.sub.,1, i.sub.t.sbsb.P.sub.,2 and for some real numbers .theta..sub.t.sbsb.P.sub.,1, .theta..sub.t.sbsb.P.sub.,2 from the interval [0,1]. If the estimated previous digital image is defined at the point (x,y), then the point (x+(t.sub.P -t.sub.k)u(x,y,t.sub.k),y+(t.sub.P -t.sub.k)v(x,y,t.sub.k)) of the gird .OMEGA..sub.t.sbsb.P.sub.,f (u(t.sub.k),v(t.sub.k)) is surrounded by the following points of the grid .OMEGA..sub.t.sbsb.P.sub.,f :

and the value I.sub.t.sbsb.P.sub.,b (x+(t.sub.P -t.sub.k)u(x,y,t.sub.k),y+(t.sub.P -t.sub.k)v(x,y,t.sub.k)) of the estimated previous digital image band is given as the following bilinear interpolation ##EQU24## of the value of the previous digital image band at the image points defined by the relations (9-8), (9-9), (9-10), (9-11).

If the estimated next digital image and the estimated previous digital image are both defined at the image point (x,y), then the value I.sub.t.sbsb.k.sub.,b (x,y,u(x,y,t.sub.k),v(x,y,t.sub.k)) of the missing digital image band at the image point (x,y) is given by the relation ##EQU25## where the value I.sub.t.sbsb.Q.sub.,b (x+(t.sub.Q -t.sub.k)u(x,y,t.sub.k),y+(t.sub.Q -t.sub.k)v(x,y,t.sub.k)) of the estimated next digital image band is given as the bilinear interpolation by the relation (9-6), while the value I.sub.t.sbsb.P.sub.,b (x+(t.sub.P -t.sub.k)u(x,y,t.sub.k),y+(t.sub.P -t.sub.k)v(x,y,t.sub.k)) of the estimated previous digital image band is given as the bilinear interpolation by the relation (9-12). If only the estimated next digital image is defined at the image point (x,y), then the value I.sub.t.sbsb.k.sub.,b (x,y,u(x,y,t.sub.k),v(x,y,t.sub.k)) of the missing digital image band at the image point (x,y) is given by the relation

where the value I.sub.t.sbsb.Q.sub.,b (x+(t.sub.Q -t.sub.k)u(x,y,t.sub.k),y+(t.sub.Q -t.sub.k)v(x,y,t.sub.k)) of the estimated next digital image band is given as the bilinear interpolation by the relation (9-6). If only the estimated previous digital image is defined at the image point (x,y), then the value I.sub.k,b (x,y,u(x,y,t.sub.k),v(x,y,t.sub.k)) of the missing digital image band at the image point (x,y) is given by the relation

where the value I.sub.t.sbsb.P.sub.,b (x+(t.sub.P -t.sub.k)u(x,y,t.sub.k),y+(t.sub.P -t.sub.k)v(x,y,t.sub.k)) of the estimated previous digital image band is given as the bilinear interpolation by the relation (9-12).

The method of computing the estimate of the velocity vector field (u(t.sub.k),v(t.sub.k)) corresponding to the missing image at the time moment t.sub.k, for each k=1, . . . K, can be described as the following multi-level resolution process:

1. For each band b=1, . . . , B, read the image band I.sub.t.sbsb.P.sub.,b into the array representing the grid .OMEGA..sub.t.sbsb.P. For each resolution level .sigma.=.sigma..sub.1, . . . , .sigma..sub.M, form the partitioning the grid .OMEGA..sub.t.sbsb.P.sup..sigma. into disjoint regions .OMEGA..sub.t.sbsb.P.sub.,f.sup..sigma., f=1, . . . , F, and then for each index g.epsilon.G representing the band b form the following: the previous filtered image function g.sup..sigma. (x,y,t.sub.P) according to the relation (4-1), the x-spatial partial derivative g.sub.x.sup..sigma. (x,y,t.sub.P) of the previous filtered image function according to the relation (4-2), and the y-spatial partial derivative g.sub.y.sup..sigma. (x,y,t.sub.P) of the previous filtered image function according to the relation (4-3).

2. For each band b=1, . . . , B, read the image band I.sub.t.sbsb.Q.sub.,b into the array representing the grid .OMEGA..sub.t.sbsb.Q. For each resolution level .sigma.=.sigma..sub.1, . . . , .sigma..sub.M, form the partitioning the grid .OMEGA..sub.t.sbsb.Q.sup..sigma. into disjoint regions .OMEGA..sub.t.sbsb.Q.sub.,f.sup..sigma., f=1, . . . , F, and then for each index g.epsilon.G representing the band b form the following: the next filtered image function g.sup..sigma. (x,y,t.sub.Q) according to the relation (4-1), the x-spatial partial derivative g.sub.x.sup..sigma. (x,y,t.sub.Q) of the next filtered image function according to the relation (4-2), and the y-spatial derivative g.sub.y.sup..sigma. (x,y,t.sub.Q) of the next filtered image function according to the relation (4-3).

3. Start with the label value F taken as the current label value f.

4. Start with the coarsest resolution level .sigma..sub.M taken as the current resolution level .sigma..

5. In the case of the current resolution level .sigma. being the coarsest resolution level .sigma..sub.M, set the initial estimate (u.sub.f,0.sup..sigma. (t.sub.k),v.sub.f,0.sup..sigma. (t.sub.k)) of the velocity vector field to some given (for example, constant) value, for every k=1, . . . , K; otherwise project the final estimate (u.sub.f,1.sup..delta.(.sigma.) (t.sub.k),v.sub.f,1.sup..delta.(.sigma.) (t.sub.k)) of the velocity vector field at the previous coarser resolution level .delta.(.sigma.)=.sigma..sub.m+1 to the current resolution level .sigma.=.sigma..sub.m to obtain the initial estimate (u.sub.f,0.sup..sigma. (t.sub.k),v.sub.f,0.sup..sigma. (t.sub.k)) of the velocity vector field at the current resolution level .sigma., for every k=1, . . . , K. Take the initial estimate (u.sub.f,0.sup..sigma. (t.sub.k),v.sub.f,0.sup..sigma. (t.sub.k)) of the velocity vector field at the current resolution level .sigma. as the current estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field at the current resolution level .sigma., for every k=1, . . . , K.

6. Form the values of the expression .parallel.(.gradient.u.sub.f.sup..sigma. (x,y,t.sub.k),.gradient.v.sub.f.sup..sigma. (x,y,t.sub.k)).parallel..sup.2, for every (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma., and for every k=1, . . . , K according to the relation (6-2).

7. For every k=1, . . . , K, the values of the following functions are formed: the optical flow function g.sub.t.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) according to the relation (5-11), for every (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma., and for every g.epsilon.G; the spatial partial with respect to the component u.sub.f.sup..sigma. (x,y,t.sub.k) derivative g.sub.tu.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the optical flow function according to the relation (5-14), for every (x,y,).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma., and for every g.epsilon.G; the spatial partial with respect to the component v.sub.f.sup..sigma. (x,y,t.sub.k) derivative g.sub.tv.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k)) of the optical flow function according to the relation (5-17), for every (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma., and for every g.epsilon.G; the expression b.sup.2 (s,.gradient.'g.sub.t.sup..sigma. (f,x,y,t.sub.k,u.sub.f.sup..sigma. (x,y,t.sub.k),v.sub.f.sup..sigma. (x,y,t.sub.k))).sup.2 according to the relation (6-3), for every (x,y).epsilon..OMEGA..sub.t.sbsb.k.sup..sigma., and for every s.epsilon.S.

8. Form the system of linear equations (6-5), for every k=1, . . . , K.

9. Obtain the improvement (.DELTA.u.sub.f.sup..sigma. (t.sub.k),.DELTA.v.sub.f.sup..sigma. (t.sub.k)) of the current estimate of the velocity vector field (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)), for every k=1, . . . , K, by iteratively solving the system of linear equations (6-5) as described in Section VII.

10. Take the improved estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k))+(.DELTA.u.sub.f.sup..sigma. (t.sub.k),.DELTA.v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field as the new current estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field, for every k=1, . . . , K.

11. Repeat step 6 to step 10 until the current estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field is improved to the desired degree, for every k=1, . . . , K. Take the resulted improved current estimate (u.sub.f.sup..sigma. (t.sub.k),v.sub.f.sup..sigma. (t.sub.k)) of the velocity vector field as the final estimate (u.sub.f,1.sup..sigma. (t.sub.k),v.sub.f,1.sup..sigma. (t.sub.k)) of the velocity vector field at the current resolution level .sigma., for every k=1, . . . , K.

12. Take the next finer resolution level .sigma..sub.m-1 as the current resolution level .sigma..

13. Repeat step 5 to step 12 until the finest resolution level .sigma.=.sigma..sub.1 is reached.

14. Take the next smaller label value f-1 as the current label value f.

15. Repeat step 4 to step 14 until the smallest label value f=1 is reached.

16. Bilinearly interpolate the final estimate (u.sub.f,1.sup..sigma. (t.sub.k),v.sub.f,1.sup..sigma. (t.sub.k)) of the finest resolution level .sigma.=.sigma..sub.1 to the grid of points .OMEGA..sub.t.sbsb.k on the image plane to obtain the estimate (u.sub.f (t.sub.k),v.sub.f (t.sub.k)) of the velocity vector field corresponding to the label value f and to the missing image at the time moment t.sub.k, for every f=1, . . . , F, and every k=1, . . . , K.

17. Combine the estimate (u.sub.f (t.sub.k),v.sub.f (t.sub.k)), f=1, . . . , F of the velocity vector field to obtain the estimate (u(t.sub.k),v(t.sub.k)) of the velocity vector field, for every k=1, . . . , K.